Usage
The reflimR package includes several functions, with
reflim(x) serving as the main function. It calls the other
important functions that can be arranged in three groups:
| group 1 | group 2 | group 3 |
|---|---|---|
| ri_hist | lognorm | adjust_digits |
| permissible_uncertainty | iboxplot | bowley |
| interpretation | truncated_qqplot | conf_int95 |
Group 1 includes the main functions that provide the user with the
final results of reflim() and help to interpret them:
ri_hist() creates a graphical output (histogram with
reference interval and density curves),
permissible_uncertainty() calculates the medical tolerance
limits of the results [4, 5], and interpretation() assesses
the medical significance of deviations from given target values [2].
Group 2 performs the three underlying statistical operations
modelling, truncation, and calculation, and
group 3 contains auxiliary functions for miscellaneous tasks. Details of
each function are available in the respective help files, which can be
addressed with the help() function or a question mark
followed by the function name (e.g. help(reflim) or
?lognorm).
The reflim function can be executed by simply calling
reflim(x), where x is a vector of numeric data to be
analyzed. For example, if the reference interval for bilirubin is to be
estimated from the values included in the livertests dataset,
the following two lines of code will suffice:
#> $stats
#> meanlog sdlog n.total n.trunc
#> 1.944 0.458 612.000 492.000
#>
#> $lognormal
#> [1] TRUE
#>
#> $limits
#> lower.lim upper.lim lower.lim.low lower.lim.upp upper.lim.low
#> 2.85 17.17 2.52 3.18 15.83
#> upper.lim.upp
#> 18.51
#>
#> $targets
#> lower.lim upper.lim lower.lim.low lower.lim.up upper.lim.low
#> NA NA NA NA NA
#> upper.lim.upp
#> NA
#>
#> $perc.norm
#> [1] 84.6
#>
#> $confidence.int
#> lower.lim.low lower.lim.upp upper.lim.low upper.lim.upp
#> 2.57 3.39 14.43 19.07
#>
#> $interpretation
#> lower.limit upper.limit
#> NA NA
#>
#> $remarks
#> [1] NAPlease note that any filtering and cleaning (e. g. for removal of
non-numeric values or duplicate results of a single individual) or
partitioning of the data (for sex, age, etc) is not included in the
reflimR package. These steps can easily be performed separately
with base R functions such as as.numeric(),
unique() or subset() to achieve the desired
results.
As shown above, the main result of the reflim function is an illustrative graphic showing a histogram of the original data with a dotted overall density curve. The estimated reference limits are displayed as dashed vertical lines and the respective medical tolerance limits as gray bars. The blue solid density curve represents the theoretical distribution of the assumed reference population, and the density curves of potential pathological outliers are shown in red.
Below the graphic, there is a list of statistical results such as the parameters of the theoretical distribution (stats), the assumed distribution model (lognormal = TRUE or FALSE), and the calculated reference limits with tolerance intervals (limits). On the bottom of the list is an estimate of the 95 percent confidence intervals for the lower and upper reference limits (confidence.int). The remaining results are set to NA if no target values were specified.
The reflim function comes with additional arguments that
default to meaningful values for a quick overview. They are accessible
via the help file (?reflim) and can be changed if needed.
For example, if you want the graphic to be nicely labeled and the
targets from the targetvalues dataset to be verified, you can
specify the main, xlab, and targets
arguments. And if you want to limit the output of the function to some
essentials such as just the interpretation of the graphic, you can
specify what you want to see with a $ after the closing parenthesis.
x <- livertests$PROT
reflim(x, main = "total protein", xlab = "g/L", targets = targetvalues[8, 3 : 4])$interpretation
#> lower.limit upper.limit
#> "within tolerance" "slightly decreased"Interpretation with traffic light colors
If target values have been specified, the other two functions of group 1 are used to check whether the calculated reference limits are within the permitted tolerance limits and assigns traffic light colors to visually indicate any deviations. Green means that the calculated limit is within the tolerance interval of the target, yellow means that the calculated limit falls outside but the two tolerance intervals overlap, and red means that there is no overlap between the tolerance intervals.
ri_hist(), permissible_uncertainty() and
interpretation() are called by the reflim
function, but can also be used independently for the evaluation of
reference limits obtained by other methods (e.g., from the
refineR package).
Three steps to the result
As already mentioned, reflimR performs an indirect three-step approach. The underlying functions form the group 2 in the above table. Here is an overview with some practical examples.
Step 1: To define an appropriate distribution model, Bowley’s quartile skewness coefficient [6] is calculated for the original and the logarithmized data, using the bowley function from group 3. The underlying algorithm is quite robust against pathological outliers, because it is calculated from the central 50 percent of the data (the “box” of Tukey’s boxplot), where the influence of pathological values is minimized [7]. If the skewness coefficient is close to 0, the distribution is symmetric and can be approximated by a normal distribution (at least in the range of interest between the 2.5th and 97.5th percentiles). If, however, the skewness coefficient is clearly positive, the distribution is right-skewed. If the shape of the curve becomes more symmetric after log-transformation of the data, the lognorm function indicates that the distribution can be better approximated by a lognormal model.
Here is an example of a symmetric distribution from the albumin data using the lognorm function:
#> $lognormal
#> [1] FALSE
#>
#> $BowleySkewness
#> normal lognormal delta
#> 0.019 -0.019 0.038The skewness coefficient of the original values is 0.019 and the shape of the black curve is not much changed after log-transformation (blue). The difference of the two skewness coefficients is 0.038, i. e. below the threshold of 0.05, which has been defined as a fixed value in the lognorm function [7].
In contrast, alanine aminotransferase is an example for an asymmetric distribution:
#> $lognormal
#> [1] TRUE
#>
#> $BowleySkewness
#> normal lognormal delta
#> 0.207 0.036 0.171Here, the original distribution curve (black) is clearly right-skewed with a skewness coefficient of 0.207, and becomes more symmetrical after logarithmization of the data (blue). The delta of the skewness coefficients is 0.171 and thus above the threshold of 0.05.
If a lognormal distribution is recommended in step 1, the values are logarithmized by the reflim function so that the normal distribution assumption can be used in the next two steps.
Step 2: Subsequently, the data set is truncated iteratively in multiple steps to remove pathological values while preserving the central 95% of the presumably non-pathological values as accurately as possible [8]. Instead of immediately truncating at the theoretical 2.5th and 97.5th percentiles (as implemented in the original package version 1.0.6), a more cautious, step-by-step truncation process has been introduced in the updated package version 1.1.0. This is done by successively adjusting the truncation threshold to the 99.5th, 99th, 98.5th, 0.98th, and 97.5th percentiles. After the last truncation step, the algorithm is continued with the final threshold of 97.5 percent until no more outliers can be removed.
Here is an application of the iboxplot function using the bilirubin dataset as an example:
The upper boxplot (white) shows a considerable amount of dots outside the whiskers, which are effectively removed with seven cycles of the iboxplot algorithm (blue boxplot). The eighth cycles confirms that the algorithm has converged:
trunc.bil$progress
#> n min max
#> cycle0 612 0.8 209.0
#> cycle1 581 2.2 24.9
#> cycle2 550 2.6 19.7
#> cycle3 541 2.7 18.6
#> cycle4 529 2.8 17.7
#> cycle5 521 2.9 17.3
#> cycle6 510 3.0 16.4
#> cycle7 492 3.1 15.3
#> cycle8 492 3.1 15.3The iboxplot function derives the estimated percentage of “normal values” from the length of the truncated vector trunc.bil$trunc, divided by the proportion of assumed non-pathological values (e.g. 0.95 as a default) :
Notably, the lower and upper truncation points are rough estimates of the reference limits (i.e. the 2.5th and 97.5th quantile of the non-pathological values):
Step 3: To improve the robustness of this estimate
even further, a normal quantile-quantile plot (Q-Q plot) is created from
the truncated data [9]. The function truncated_qqplot()
determines the mean and standard deviation from the intercept and the
slope of the regression line and extrapolates the 2.5th and 97.5th
percentiles from the equation mean ± 1.96 sd. These percentiles
represent the estimated reference interval, and thus the final
result.
#> $result
#> meanlog sdlog lower.lim upper.lim
#> 1.944 0.458 2.850 17.170
#>
#> $lognormal
#> [1] TRUEIf the calculation was carried out with logarithmized values, the resulting reference limits are automatically delogarithmized.
As already mentioned in step 2, the truncation points of step 2 are close to the final result of step 3, but they are not identical. The reason is that extrapolating the respective percentiles from the linear Q-Q plot obviously stabilizes the result and makes it more robust against minor overlaps with pathological values at the ends of the regression line (see above). To achieve this, only the central 90 percent of the Q-Q plot are used for the calculation of the regression line.
