This report is automatically generated with the R package knitr (version 1.40) by .

library(fitdistrplus) # fits data to distributions - fitdistcens() & fitdist() - load first so it doesn't mask dplyr::select
library(tidyverse)    # functions to work with tibbles


# FRACTION SUBSTITUTION OF MDL or RL VALUE WHEN RESULT IS BLANK & RESULT QUAL CODE IS DNQ OR ND #
substituteNDDNQ <- function(dataframe, QualCodeColumn, MDL.RL.fraction=.5){ # QualCodeColumn <- quo(NewQualCode) - to test funtion line-wise
  qualColumn <- enquo(QualCodeColumn)

  randResult <- dataframe %>%
    mutate(Censored   = case_when(!!qualColumn=='ND'                  ~ T,
                                  !!qualColumn=='DNQ' & is.na(Result) ~ T,            # Label DNQ as censored only if Result is blank
                                   TRUE ~ F),
           Result = case_when(!!qualColumn=='ND' & !is.na(MDL)    ~ map_dbl(.x=MDL, ~ .x * MDL.RL.fraction),             # '~' denotes  to use a formula, '.x' is the vector - in this case, the MDL column
                              !!qualColumn=='ND'                  ~ map_dbl(.x=RL,  ~ .x * MDL.RL.fraction),
                              !!qualColumn=='DNQ' & is.na(Result) ~ map2_dbl(.x=MDL, .y=RL, ~mean(c(.x, .y))),
                              TRUE ~ Result)
    )
  return(randResult)
}


# RANDOMIZE RESULT VALUE WHEN RESULT IS BLANK & RESULT QUAL CODE IS DNQ OR ND #
randomNDDNQ <- function(dataframe, QualCodeColumn){ # QualCodeColumn <- quo(NewQualCode) - to test funtion line-wise
  qualColumn <- enquo(QualCodeColumn)

  set.seed(124)
  suppressWarnings(
    randResult <- dataframe %>%
      mutate(
        # Label ND/DNQs as being censored
        Censored   = case_when(!!qualColumn=='ND'                  ~ T,
                               !!qualColumn=='DNQ' & is.na(Result) ~ T,            # Label DNQ as censored only if Result is blank
                               TRUE ~ F),
        Result = case_when(!!qualColumn=='ND' & !is.na(MDL)    ~ map_dbl(.x=MDL, ~runif(1, min=0, max=.x)),
                           !!qualColumn=='ND'                  ~ map_dbl(.x=RL , ~runif(1, min=0, max=.x)),
                           !!qualColumn=='DNQ' & is.na(Result) ~ map2_dbl(.x=MDL, .y=RL, ~runif(1, min=.x, max=.y)),
                           TRUE ~ Result)
             )
  )
  return(randResult)
}


# ROBUST REGRESSION ORDER OF STATISTICS (ROS) ON NON DETECTS (MAYBE BEST ALTERNATIVE TO RANDOMIZING RESULT OR 1/2 MDL SUBSTITUTION)
# ROS using NADA
library(NADA)
# https://www.itrcweb.org/gsmc-1/Content/GW%20Stats/5%20Methods%20in%20indiv%20Topics/5%207%20Nondetects.htm
# https://www.itrcweb.org/gsmc-1/Content/Resources/Unified_Guidance_2009.pdf#15.4
nadaNDDNQ <- function(dataframe, QualCodeColumn, forwardT="log", reverseT="exp"){ # QualCodeColumn <- quo(NewQualCode) - to test funtion line-wise
  qualColumn <- enquo(QualCodeColumn)

  rosResult <- dataframe %>%
    mutate(Censored = case_when(!!qualColumn=='ND'                  ~ T,
                                !!qualColumn=='DNQ' & is.na(Result) ~ T,
                                TRUE ~ F),
           Result = case_when(!!qualColumn=='ND' & !is.na(MDL)    ~ MDL,
                              !!qualColumn=='ND'                  ~ RL,
                              !!qualColumn=='DNQ' & is.na(Result) ~ RL,
                              TRUE ~ Result)
    ) %>%
    arrange(Result) # so modeled results are included in same order as original values
  nadaros <- NADA::ros(rosResult %>% pull(Result), rosResult %>% pull(Censored), forwardT=forwardT, reverseT=reverseT)
  nadaros <- as.data.frame(nadaros)

  if(isFALSE( identical(rosResult$Result, nadaros$obs) )){stop('ROS modeled results do not line up with the Observed Results')}

  #Replace original values with modeled results 
  rosResult <- rosResult %>%
    mutate(Result = nadaros$modeled)

  return(rosResult)
}


# MAXIMUM LIKLIHOOD ESTIMATION (MLE) ON NON DETECTS USING fitdistrplus::fitdistcens & fitdistrplus::fitdist iterations to improve AIC model value
# https://cran.r-project.org/web/packages/fitdistrplus/vignettes/paper2JSS.pdf
# library(fitdistrplus) # fits data to distributions - fitdistcens() & fitdist() - needs to be run before library(tideyverse) so it doesn't mask dplyr::select
# library(extraDistr)  # provides truncated normal distribution (tnorm) model - stopped using due to tnorm() errors in NDDNQ estimate selection in setNDDNQtoDistFit() 
library(actuar)       # provides additional distribution models like llogis, pareto, burr - load after library(extraDistr) so actuar::pareto functions aren't masked

# This is my approach of testing multiple distributions ('burr', 'gamma', 'llogis', 'lnorm', 'pareto', 'weibull') to fit data, then fill in ND & DNQ based on distribution model
# fitdistcens() makes a first guess at the distribution coefficients, then I use fitdist() in 2 more iterations to refine coefficients (to save processing time, the 2nd iteration is only performed if the AIC improved in the 1st iteration)
# Overall, my approach outperforms 1/2 Substitution, Random Value, & NADA ROS in Monte Carlo testing but I am challenged with selecting the best fitted distribution
# In attempts to select best distribution model, I tested choosing the model with the best RMSE, MAPE, and Correlation (all comparing the random selected modeled results with the distribution model results), and AIC. No parameter consistently selected the best distribution model.
# I chose AIC as the model selection parameter because it is calculated by fitdist() and is consistently recommended in the literature.
# Because a parameter did not consistently select the best distribution model, The Monte Carlo test evaluated which models tend to have the best fit overall.
# The Monte Carlo test results show llogis, then lnorm, then weibull perform best. Order of returned distribution results default to llogis, then lnorm, then weibull, then best to worst AIC value.
# User can select first result by default with some confidence that the best modeled ND & DNQ values were chosen. Or, they can look at the other distribution results and manually choose which results to use.

# SECONDARY FUNCTIONS TO SUPPORT MAIN FUNCTION robinNDDNQ()#
# Opposite of %in%
`%are not%` <- function (x, table) is.na(match(x, table, nomatch=NA_integer_)) #https://stackoverflow.com/questions/5831794/opposite-of-in

# Collapses vector to be like a grammatical list (i.e., seperated by commas, with individual quotes, and last element "and")  #
# Needs to be wrapped in a cat() to print properly; e.g., cat(grammaticList(c('a','b', 'c')))
grammaticList <- function(vec) sub(",\\s+([^,]+)$", ", and \\1", ifelse(length(vec) == 2, paste(vec, collapse=" and "), paste(vec, collapse=", ")))

# Replace ND & DNQ values with values from fitted distribution using probability of <MDL & <RL values
setNDDNQtoDistFit <- function(quals, distname, coef.est, fit, dataframe_cens, dataframe_mle){
  # QualCodeColumn <- quo(NewQualCode) # to test function line-wise
  # qualColumn <- enquo(QualCodeColumn)

  # To generate possible concentrations Summarize quals so DNQs with MDL=RL are managed as NDs
  # (prevents instances of random selection having no concentrations to choose from)
  # qual_summary <- quals %>% 
  #   mutate(!!qualColumn := if_else(RL == tempMDL, 'ND', !!qualColumn)) %>% # change DNQs with RL=MDL to NDs
  #   mutate(RL = if_else(!!qualColumn == 'ND', tempMDL, RL)) %>%            # for NDs, replace -88 in RL column with MDL value
  #   group_by(!!qualColumn, tempMDL, RL) %>%                                # group similar MDL/RL values
  #   summarise(n=sum(n), .groups='drop') %>%                                # sum the total n * use this to determine the number estimated concentrations that are needed fo random selection
  #   arrange(desc(n))

  # Randomly select distribution model concentrations
  for (multiplier in 2:10) { # initially set multiplier low to generate # of possible concentrations, if random sample() generates error due to not enough concentrations, try a larger multiplier next  
    fitConcs <- NULL
    # generate distribution estimated Concentrations less than MDL or RL, or between MDL & RL based on fitted distribution
    for (row in 1:nrow(quals)) {
      if (quals[['tempQual']][row] == 'ND') { # ND or DNQ without different RL/MDL
        prob.string.MDL <- paste0('p', distname,'(', quals$tempMDL[row], ', ', paste(names(coef.est) ,'=', coef.est, collapse=', '), ')') # create dynamic chr string to evaluate probability of MDL
        prob.MDL <- eval(parse(text=prob.string.MDL))                                                                                            # get probablity for fitted Results less than MDL
        subfitConcs <- quantile(fit, seq(0, prob.MDL, length=quals$n[row] * multiplier))                                                   # requires library(fitdistrplus) to work properly; for help ?fitdistrplus::quantiles; returns list of quantiles (i.e. possible concentrations) for specified probabilities
        fitConcs <- rbind(fitConcs, data.frame(quantile=unlist(subfitConcs$quantiles), prob=subfitConcs$probs))

      } else { # DNQ and need to select between MDL and RL
        prob.string.MDL <- paste0('p', distname,'(', quals$tempMDL[row], ', ', paste(names(coef.est) ,'=', coef.est, collapse=', '), ')') #create dynamic chr string to evaluate probability of MDL
        prob.MDL <- eval(parse(text=prob.string.MDL))                                                                                            # get probablity for fitted Results less than MDL
        prob.string.RL <- paste0('p', distname,'(', quals$tempRL[row], ', ', paste(names(coef.est) ,'=', coef.est, collapse=', '), ')')       # create dynamic chr string to evaluate probability of RL
        prob.RL <- eval(parse(text=prob.string.RL))                                                                                              # get probablity for fitted Results less than RL
        subfitConcs <- quantile(fit, seq(prob.MDL, prob.RL, length=quals$n[row] * multiplier))                                             # requires library(fitdistrplus) to work properly; for help ?fitdistrplus::quantiles; returns list of quantiles (i.e. possible concentrations) for specified probabilities
        fitConcs <- rbind(fitConcs, data.frame(quantile=unlist(subfitConcs$quantiles), prob=subfitConcs$probs))
      }
    }
    fitConcs <- fitConcs %>% distinct(quantile, .keep_all=T) %>% arrange(quantile)

    # Randomly select NDDNQ estimates
    for (row in 1:nrow(quals)) {
      if (quals[['tempQual']][row] == 'ND') {
        # Randomly sample from distribution modeled concentrations as estimate for NDs
        quantiles <- fitConcs %>%
          filter(quantile > 0 & signif(quantile, 7) < quals$tempMDL[row])
        # observed cases where fitConcs(MDL value) != MDL, so fitConcs(MDL value) > MDL was T but should be F, so included signif() as safety measure to exclude MDL & RL value from random sample options

        set.seed(124)                                                         # set seed so same concentrations are selected if sample is run again (i.e., so results are reproducible)
        capture.output( mle_prob_sample <- try(sample(x=quantiles$quantile, size=quals$n[row], replace=F),  # sample fitted concentrations less than MDL without replacement (if "prob=quantiles$prob" is added to base selection on probabilities, Monte Carlo tests show much worse performance)
                                               silent=T), file='NUL' )                             # use capture.output to stop printing cat warning
        if(inherits(mle_prob_sample, "try-error")) next                # if sample() results in an error, try next multiplier to generate possible concentrations

      } else { # qualColumn == 'DNQ'
        # Randomly sample from distribution modeled concentrations as estimate for DNQs
        if (quals$tempMDL[row] == quals$tempRL[row]) {
          # Sample Concentrations less than RL
          quantiles <- fitConcs %>%
            filter(quantile > 0 & signif(quantile, 7) < quals$tempRL[row])
          # observed cases where fitConcs(MDL value) != MDL, so fitConcs(MDL value) > MDL was T but should be F, so included signif() as safety measure to exclude MDL & RL value from random sample options


          set.seed(124)                                                         # set seed so same concentrations are selected if sample is run again (i.e., so results are reproducible)
          capture.output( mle_prob_sample <- try(sample(x=quantiles$quantile, size=quals$n[row], replace=F),  # sample fitted concentrations less than MDL without replacement (if "prob=quantiles$prob" is added to base selection on probabilities, Monte Carlo tests show much worse performance)
                                                 silent=T), file='NUL' )                             # use capture.output to stop printing cat warning
          if(inherits(mle_prob_sample, "try-error")) next                # if sample() results in an error, try next multiplier to generate possible concentrations

        } else {
          # Sample Concentrations between MDL & RL
          quantiles <- fitConcs %>%
            filter(signif(quantile, 7) > quals$tempMDL[row] & signif(quantile, 7) < quals$tempRL[row])
          # observed cases where fitConcs(MDL value) != MDL, so fitConcs(MDL value) > MDL was T but should be F, so included signif() as safety measure to exclude MDL & RL value from random sample options

          set.seed(124)                                                         # set seed so same concentrations are selected if sample is run again (i.e., so results are reproducible)
          capture.output( mle_prob_sample <- try(sample(x=quantiles$quantile, size=quals$n[row], replace=F, prob=quantiles$prob),  # sample fitted concentrations less than MDL without replacement (if "probs=probs" is added to base selection on probabilities, Monte Carlo tests show much worse performance)
                                                 silent=T), file='NUL' )                             # use capture.output to stop printing cat warning
          if(inherits(mle_prob_sample, "try-error")) next                # if sample() results in an error, try next multiplier to generate possible concentrations
        }
      }
      # Remove selected sample from fitConcs pool for next row of quals random selection
      fitConcs <- fitConcs %>%
        anti_join(., data.frame(quantile=mle_prob_sample), by='quantile')

      # Get indices of ND/DNQ instances
      indices <- dataframe_cens %>%
        mutate(row.id = row_number()) %>%
        filter(tempQual == quals[['tempQual']][row],
               tempMDL  == quals$tempMDL[row],
               tempRL   == quals$tempRL[row],
               Censored == T) %>%
        pull(row.id)

      # Replace ND/DNQ occurrences with randomly selected modeled concentrations
      dataframe_mle$Result <- replace(dataframe_mle$Result, indices, mle_prob_sample)
    }
  }
  return(dataframe_mle)
}

# MAIN FUNCTION #
# My function to estimate ND & DNQs based on fitted distribution model
robinNDDNQ <- function(dataframe, QualCodeColumn){
  # QualCodeColumn <- quo(ResultQualCode) # to test function line-wise
  qualColumn <- enquo(QualCodeColumn)

  # Setup dataframe with censored, left, and right columns for fitdistcens()
  dataframe_cens <- dataframe %>%
    mutate(origResult = Result,                                                       # Keep track of original result column for future reference
           # Create temp Qual, RL, & MDL columns for easier coding later 
           tempMDL   = if_else(!!qualColumn == 'ND' & is.na(MDL), RL, MDL), # For NDs & MDL value is blank, enter RL value in MDL column (RL substutiting for MDL in these cases)
           tempQual  = if_else(!!qualColumn == 'DNQ' & RL == tempMDL, 'ND', !!qualColumn),          # change DNQs to NDs when RL=MDL
           tempRL    = if_else(tempQual == 'ND', tempMDL, RL),              # For NDs, enter RL value into MDL column (helps deal with missing RL values)
           # Label ND/DNQs as being censored
           Censored   = case_when(tempQual=='ND'                  ~ T,
                                  tempQual=='DNQ' & is.na(Result) ~ T,            # Label DNQ as censored only if Result is blank
                                  TRUE ~ F),
           # Replace censored Results with MDL/RL values (also sets up 'right' censor (max value))
           Result     = case_when(tempQual=='ND' & !is.na(MDL)    ~ MDL,
                                  tempQual=='ND'                  ~ RL,
                                  tempQual=='DNQ' & is.na(Result) ~ RL,           # Right censor DNQ (concentration between MDL & RL) only if Result is blank
                                  TRUE ~ Result),
           # Set up left censor column (min value), use NA_real for NDs since they have no lower limit
           left       = case_when(tempQual=='ND'                      ~ NA_real_,
                                  tempQual=='DNQ' & is.na(origResult) ~ MDL,      # Left censor DNQ (concentration between MDL & RL) only if Result is blank
                                  TRUE ~ Result),
           # Set up right censor column using the work done for the Result column
           right      = Result) %>%
    as.data.frame()
  dataframe_mle <- dataframe_cens %>%
    select(-left, -right) %>%    #removes left & right columns from output
    as.data.frame()

  # Summarize & count different types of Censored conditions
  dnqs <- dataframe_cens %>% filter(tempQual == 'DNQ' & is.na(origResult)) %>% group_by(tempQual, tempMDL, tempRL) %>% summarise(n=n(), .groups='drop') %>%  filter(!is.na(tempMDL) & !is.na(tempRL))
  nds  <- dataframe_cens %>% filter(tempQual == 'ND')                      %>% group_by(tempQual, tempMDL, tempRL) %>% summarise(n=n(), .groups='drop') %>%  filter(!is.na(tempMDL) & !is.na(tempRL))
  quals <- rbind(nds, dnqs) %>%
    group_by(tempQual, tempMDL, tempRL) %>%  # group similar MDL/RL values
    summarise(n=sum(n), .groups='drop') %>%  # sum the total n * use this to determine the number estimated concentrations that are needed fo random selection
    arrange(desc(n))

  # Estimate population mean & sd using sample values to set initial estimate of distribution coefficients
  estNDDNQ   <- randomNDDNQ(dataframe_cens %>% mutate(Result=origResult), tempQual) %>% pull(Result)
  estMean    <- mean(estNDDNQ)
  estLogMean <- mean(log(estNDDNQ))
  estSD      <- sd(estNDDNQ)
  estLogSD   <- sd(log(estNDDNQ))

  # Only continue if NDs/DNQs are present in the sample set 
  if (nrow(quals) > 0) {

    # Identify distributions to test and initial estimate of distribution coefficients
    distNames <- c('burr', 'gamma', 'llogis', 'lnorm', 'pareto', 'weibull')  # stopped using tnorm due to errors during estimated NDDNQ selection in in setNDDNQtoDistFit()
    distStart <- list(burr    = list(shape1=0.3, shape2=1, rate=1),
                      gamma   = list(shape=estMean^2/estSD^2, scale=estSD^2/estMean),
                      llogis  = NULL,
                      lnorm   = list(meanlog=estLogMean, sdlog=estLogSD),
                      tnorm   = list(mean=estMean, sd=estSD),             # left tnorm here just in case I decide to start using it again
                      pareto  = list(shape=1, scale=500),
                      weibull = list(shape=(estSD/estMean)^-1.086, scale=estMean/gamma(1+1/(estSD/estMean)^-1.086)))
    fixArg <- list(burr=NULL, gamma=NULL, llogis=NULL, lnorm=NULL, tnorm=list(a=0, b=Inf), pareto=NULL, weibull=NULL) # left tnorm here just in case I decide to start using it again

    # Fit the sample set to each distribtuion                         
    mleResults <- list()
    for (distname in distNames) {
      #FIRST - fit distribution to censored data using Maximum Likelihood Estimation (MLE)
      capture.output( first.fit <- try(fitdistcens(dataframe_cens, distname, fix.arg=fixArg[[distname]], start=distStart[[distname]], lower=c(0, 0)),
                                       silent=T), file='NUL' ) # use capture.output to stop printing cat warning 
      if(inherits(first.fit, "try-error")) next                # if fitdistcens() results in an error, goes to next distribution

      # Use custom function setNDDNQtoDistFit() to replace ND/DNQs with randomly selected estimates from distribution fit 
      dataframe_mle1 <- setNDDNQtoDistFit(quals, distname, coef.est=coef(first.fit), fit=first.fit, dataframe_cens, dataframe_mle)


      #SECOND - fit distribution to modeled censored data & refine coefficient estimates
      capture.output( second.fit <- try(fitdistrplus::fitdist(dataframe_mle1$Result, distname, fix.arg=fixArg[[distname]], start=distStart[[distname]], lower=c(0, 0)),
                                        silent=T), file='NUL' )
      if(inherits(second.fit, "try-error")) second.fit <- first.fit

      # Use custom function setNDDNQtoDistFit() to replace ND/DNQs with randomly selected estimates from second distribution fit 
      dataframe_mle2 <- setNDDNQtoDistFit(quals, distname, coef.est=coef(second.fit), fit=second.fit, dataframe_cens, dataframe_mle)


      #If first.fit has lower aic, store results, otherwise do one more iteration - lower aic is better
      if(first.fit$aic < second.fit$aic){  # store dataframe_mle1 (original dataframe with fitted results) and fit model
        dataframe_mle <- dataframe_mle1
        final.fit     <- first.fit
      } else { # if aic improved, do another iteration to see if aic improves more
        #THIRD - fit distribution to modeled censored data & refine coefficient estimates
        capture.output( third.fit <- try(fitdistrplus::fitdist(dataframe_mle2$Result, distname, fix.arg=fixArg[[distname]], start=distStart[[distname]], lower=c(0, 0)),
                                         silent=T), file='NUL' )
        if(inherits(third.fit, "try-error")) third.fit <- second.fit

        # Use custom function setNDDNQtoDistFit() to replace ND/DNQs with randomly selected estimates from third distribution fit 
        dataframe_mle3 <- setNDDNQtoDistFit(quals, distname, coef.est=coef(third.fit), fit=third.fit, dataframe_cens, dataframe_mle)

        #If second.fit has lower aic, store results, otherwise store third.fit - lower aic is better
        if(second.fit$aic < third.fit$aic){ # store dataframe_mle2 (original dataframe with fitted results) and fit model
          dataframe_mle <- dataframe_mle2
          final.fit     <- second.fit
        } else {                            # store dataframe_mle3 (original dataframe with fitted results) and fit model
          dataframe_mle <- dataframe_mle3
          final.fit     <- third.fit
        }
      }
      # STORE ORIGINAL DATAFRAME WITH MODELED CENSORED VALUES, FINAL FIT MODEL, & AIC
      mleResults[[distname]][['fitted']]   <- dataframe_mle %>% select(-tempQual, -tempRL, -tempMDL)
      mleResults[[distname]][['fitmodel']] <- final.fit
      mleResults[[distname]][['aic']]      <- final.fit$aic # Store aic for distribution model selection, lower aic is better
    }


    # ORDER RESULTS # - Ideally best fitted distribution would be listed first, but could not find metric to consistently rank best distribution fit.
    # So went with AIC, which was also used for the fitdistrplus::fitdist() iteration selection 
    mleResults <- mleResults[order(sapply(mleResults, `[[`, i = 'aic'), decreasing=F)] #order list so best (lowest) AIC is first
    dist.tested <- names(mleResults)
    # Based on Monte Carlo results, best order seems to be #1 llogis, #2 weibull, and #3 gamma (gamma and lnorm essentially tied but went with gamma as it has been shown to be a good representation of environmental data)
    # Working backward reorder list so llogis is 1st, weibull 2nd, & gamma 3rd
    if('gamma' %in% dist.tested){
      dist.tested <- c('gamma', dist.tested[dist.tested %are not% 'gamma'])     # make gamma first
      mleResults <- mleResults[dist.tested] }
    if('weibull' %in% dist.tested){
      dist.tested <- c('weibull', dist.tested[dist.tested %are not% 'weibull'])   # make weibull first, this will make gamma second
      mleResults <- mleResults[dist.tested] }
    if('llogis' %in% dist.tested){
      dist.tested <- c('llogis', dist.tested[dist.tested %are not% 'llogis']) # make llogis first, this will make weibull second, and gamma third
      mleResults <- mleResults[dist.tested] }


    cat('Tested', grammaticList(dist.tested), '\nResults for tested distributions are returned as a list. \n',
        'Order of distributions defaults to llogis, then weibull, then gamma, then best to worst AIC value.\n',
        'To retrieve dataframe with estimated ND/DNQ values use "<varname>$<distname>$fitted"    or      "<varname>[[1]][[1]]"\n',
        'To view goodness of fit graphs                use "plot(<varname>$<distname>$fitmodel)" or "plot(<varname>[[1]][[2]])"\n',
        'Unable to test', grammaticList(distNames[distNames %are not% dist.tested]), 'distributions due to fitdistcen() errors.\n')
    #View(mleResults$llogis$fitted) 
    #plot(mleResults$llogis$fitmodel)
    #summary(mleResults$llogis$fitmodel)
  } else {
    mleResults <- dataframe
    cat('No ND or DNQ data found in dataset. Returned original dataframe.')
  }

  return(mleResults)
}


# # CODE TO TEST FUNCTIONS #
# # Random Data Setup
# mean <- 3
# sd   <- 2
# 
# shape <- mean^2/sd^2
# scale <- sd^2/mean
# set.seed(124)
# origdata <- data.frame(Result=rgamma(n=20, shape=shape, scale=scale)) %>% 
#   arrange(Result)
# 
# 
# # Mock Censored Data
# rl  <- 1.5
# mdl <- 0.8
# yourdata <- origdata %>% 
#   mutate(QualCode = case_when(Result <= mdl ~ 'ND',
#                               Result <= rl ~ 'DNQ',
#                               TRUE ~ '='),
#          Result = case_when(Result <= mdl ~ NA_real_,
#                             Result <= rl ~ NA_real_,
#                             TRUE ~ Result),
#          MDL = mdl,
#          RL = rl)
# 
# ## Required columns are "Result", "MDL", "RL", & a QualCodeColumn (you can specify the QualCodeColumn name in the equation below)
# 
# 
# # Estimate NDDNQs - see "functions_estimate NDDNQ values.R" more details
# results <- robinNDDNQ(dataframe=yourdata, QualCodeColumn=QualCode)
# results$llogis$fitted         # or results[[1]]$fitted
# plot(results$llogis$fitmodel) # or plot(results[[1]]$fitmodel)
# 
# 
# ## Other functions for substituting or estimating NDDNQs in "functions_estimate NDDNQ values.R"
# substituteNDDNQ(dataframe=yourdata, QualCodeColumn=QualCode, MDL.RL.fraction=.5)
# 
# randomNDDNQ(dataframe=yourdata, QualCodeColumn=QualCode)
# 
# nadaNDDNQ(dataframe=yourdata, QualCodeColumn=QualCode)

The R session information (including the OS info, R version and all packages used):

    sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=C                            
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## system code page: 65001
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.8.0    plotly_4.10.0      readxl_1.4.1       actuar_3.3-0      
##  [5] NADA_1.6-1.1       forcats_0.5.2      stringr_1.4.1      dplyr_1.0.9       
##  [9] purrr_0.3.4        readr_2.1.2        tidyr_1.2.0        tibble_3.1.8      
## [13] ggplot2_3.3.6      tidyverse_1.3.2    fitdistrplus_1.1-8 survival_3.4-0    
## [17] MASS_7.3-58.1     
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-45     assertthat_0.2.1    digest_0.6.29       utf8_1.2.2         
##  [5] R6_2.5.1            cellranger_1.1.0    backports_1.4.1     reprex_2.0.2       
##  [9] evaluate_0.16       highr_0.9           httr_1.4.4          pillar_1.8.1       
## [13] rlang_1.0.5         lazyeval_0.2.2      googlesheets4_1.0.1 data.table_1.14.2  
## [17] rstudioapi_0.14     Matrix_1.5-1        rmarkdown_2.16      splines_4.2.2      
## [21] googledrive_2.0.0   htmlwidgets_1.5.4   munsell_0.5.0       broom_1.0.1        
## [25] compiler_4.2.2      modelr_0.1.9        xfun_0.32           pkgconfig_2.0.3    
## [29] htmltools_0.5.3     tidyselect_1.1.2    viridisLite_0.4.1   fansi_1.0.3        
## [33] crayon_1.5.1        tzdb_0.3.0          dbplyr_2.2.1        withr_2.5.0        
## [37] grid_4.2.2          jsonlite_1.8.0      gtable_0.3.1        lifecycle_1.0.1    
## [41] DBI_1.1.3           magrittr_2.0.3      scales_1.2.1        writexl_1.4.0      
## [45] cli_3.3.0           stringi_1.7.8       fs_1.5.2            xml2_1.3.3         
## [49] ellipsis_0.3.2      generics_0.1.3      vctrs_0.4.1         expint_0.1-7       
## [53] tools_4.2.2         glue_1.6.2          hms_1.1.2           fastmap_1.1.0      
## [57] yaml_2.3.5          colorspace_2.0-3    gargle_1.2.0        rvest_1.0.3        
## [61] knitr_1.40          haven_2.5.1
    Sys.time()
## [1] "2023-12-29 08:24:38 PST"