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"