*! version 1.0 19 Jan 2023

Title

ros -- Regression order statistics - Estimating upper reference
bounds for a dataset with possibly non-detectable/censored positive values and
possibly contaminated in the upper end.

Syntax

ros
varlist(max=1
numeric)
[if]
[,
options]

Optional

censor A variable indicating whether a value is censored (1) or
not (0)

noquietly Show regression output and save generated variables if
set

percentiles A numlist of values between 0 and 100. Specifies the
percentiles to estimate. The values will be sorted. Default values are 50 75 90 95 99

scatter Generates a qnorm scatter plot as model control.
The observed values should either be on a straight line.
Or at least the first part should be.

rsqrtheta Generates a line plot of adjusted R squares by thetas
(Box-Cox transformation). The thetas with the highest R squares are the best
Box-Cox transformations.

theta(#) for the choice of Box-Cox transformation, default = 1.

start(#) Start value for the theta's in the rsqrtheta graph.
The default value is 0.

step(#) Step value for the theta's in the rsqrtheta graph.
The default value is 0.25.

end(#) end value for the theta's in the rsqrtheta graph.
The default value is 2.

twoway options Add twoway graph options to the qnorm
scatter plot

Description

The command ros is for estimating upper reference bounds for a
dataset with possibly non-detectable/censored values and possibly
contaminated in the upper end.

The upper reference bounds are from the mean and standard deviation from
regressing the observed values on the empirical (normal) z-values.
The mean and the standard deviation estimated from ros are the intercept
and slope of the regression.
Box-Cox optimization is done by choosing a transformation with a high R square
from the regressions.

A dataset can have non-detectable/censored values to the left and can also be
contaminated to right by non-healthy participants, so it might be only a
small linear part of the qnorm scatter plot that is acceptable for
estimation.
This means that some of the estimated upper bounds are extrapolations from the
acceptable part.
The estimated, not extrapolated, upper bounds, are similar to the empirical
percentiles.

The Ideas behind the command ros were demonstrated at the
2022 Northern European Stata Conference on 12 October in Oslo.
The ros is inspired by the work in the references.
However, more tools for model control, the rsquare for goodness of fit,
the rsquare theta plot, and the qnorm scatter plot are available.

The only modification from the presentation is that the Box-Cox
transformation here is defined as: bct(x; theta) = x^theta / abs(theta)
if theta != 0 and bct(x; 0) = log(x).
This modification means that negative thetas can be included in the Box-Cox
transformations.

Examples

Below are some examples of how to assess a dataset using ros.
It is worth noting that when using ros there is no correct model,
but sometimes there is a useful model.

Real data example

The first dataset is from Huston (2009). It is used to show that using
the log as a default transformation is a mistake.

Using log as a default transformation is what has been recommended in, e.g.,
Helsel (2010) and Hewett and Ganser (2009)

Get the data (Can retrieved together with this package).

`". use "savona (NADA).dta", clear"'

Get a description of the dataset.

`". notes"'

Rescale the measured concentrations with a factor 1000.

`". replace concentration = 1000 * concentration"'

Use the log transformation as done in Huston (2009), Helsel (2010),
and Hewett and Ganser (2009).

`". ros concentration, censor(censored) theta(0) rsqrtheta scatter"'

Adjusted Rsquared is 0.8653

------------------------------------------------
Estimate Empirical
------------------------------------------------
theta(0.000)
P50% 2.85 2.00
P75% 4.77 6.00
P90% 7.59 7.00
P95% 10.01 8.35
P97.5% 12.74 9.00
P99% 16.85 9.00
------------------------------------------------

--------------------------------------
R2 adj theta optimal
--------------------------------------
0.865 0.000 .
0.881 0.250 .
0.894 0.500 .
0.906 0.750 .
0.915 1.000 .
0.921 1.250 .
0.923 1.500 1
0.921 1.750 .
0.916 2.000 .
--------------------------------------

There is a relatively poor fit (rsquare = 0.8653).
This is also seen in the differences between the estimated (from ros)
and the empirical (from centile).
Note that the best is with a rsquare at 0.923 at theta equal to 1.5.

Looking at the qnorm scatter plot, there is no indication of a
mixture of distributions.

The Box-Cox transformation 1.5 is tried.

`". ros concentration, censor(censored) theta(1.5) rsqrtheta scatter"'

Adjusted Rsquared is 0.9228

------------------------------------------------
Estimate Empirical
------------------------------------------------
theta(1.500)
P50% 3.42 2.00
P75% 5.61 6.00
P90% 7.26 7.00
P95% 8.16 8.35
P97.5% 8.90 9.00
P99% 9.73 9.00
------------------------------------------------

The estimated and empirical percentiles are similar, although the
empirical median and the empirical 99% percentile are underestimated.

The estimated percentiles using the log transformation appear to be too high.

Simulated data example 1

A log-normal data set of size 200 is created with a mean of 2 and a
standard deviation of 1 for the log-transformed data,
see boxcoxsim.

`". boxcoxsim, n(200) nd(0) theta(0) mean(4) sd(1) outlierpct(0) clear"'

To visualize the data set.

`". hist yc, norm ylabel(none) xtitle(normal values)"'

Run the ros command to assess the data

`". ros yc, censor(censored) rsqrtheta scatter start(-1)"'

The adjusted rsquare is low.

From the log window and the rsqrtheta graph, one can see that the log
transformation (theta = 0) is best.

The curve in the qnorm scatter plot is a typical log curve
indicating the need for the log transformation.

The log transformation is used in the ros command.

`". ros yc, censor(censored) rsqrtheta scatter start(-1) theta(0)"'

The rsqrtheta graph indicates that the log transformation (theta = 0)
is best.

The curve in the qnorm scatter plot is linear.

In the log window, there is quite a difference between the estimated and
the empirical percentiles. In this case, the estimated percentiles should be
used for upper bounds.

Simulated data example 2

A log-normal data set of size 200 is created with a mean of 1.5 and a
standard deviation of 0.3 for the power 4 transformation data,
see boxcoxsim.

`". boxcoxsim, n(200) nd(0) theta(4) mean(1.5) sd(0.3) outlierpct(0) clear"'

To visualize the data set.

`". hist yc, norm ylabel(none) xtitle(normal values)"'

Run the ros command to assess the data

`". ros yc, censor(censored) rsqrtheta scatter start(-1)"'

The curve in the qnorm scatter plot is not linear.

The rsqrtheta graph identifies the correct inverse transformation
(theta = 0.25) as optimal.

Run the ros command with the optimal theta (0.25).

`". ros yc, censor(censored) rsqrtheta scatter start(-1) theta(0.25)"'

The curve in the qnorm scatter plot is now linear.

The rsqrtheta graph identifies the correct inverse transformation
(theta = 0.25) as optimal.

Simulated data example 3

A log-normal data set of size 200 is created with a mean of 5 and a
standard deviation of 1 for the log-transformed data,
see boxcoxsim.

There are 55% of the observations below the lower detection line.

`". boxcoxsim, n(200) nd(55) theta(0) mean(4) sd(1) outlierpct(0) clear"'

To visualize the data set.

`". hist yc, norm ylabel(none) xtitle(normal values)"'

Run the ros command to assess the data

`"ros yc, censor(censored) rsqrtheta scatter start(-1)"'

The curve in the qnorm scatter plot is curved like a log function.

The rsqrtheta graph identifies the correct inverse transformation
(theta = 0, ie the log) as optimal.

Run the ros command with the optimal theta (0).

`"ros yc, censor(censored) rsqrtheta scatter theta(0) start(-1)"'

This verifies our assessments.

Simulated data example 4

A normal data set of size 200 is created with a mean of 5 and a standard
deviation of 1, see boxcoxsim.

There are 40% of the observations below the lower detection line.

The upper 30% are outliers and the outlier distribution is normal with a mean of 11
and a standard deviation of 2.

`". boxcoxsim, n(200) nd(40) mean(5) sd(1) outlierpct(30) omean(11) osd(2) clear"'

Run the ros command to assess the data

`". ros yc, censor(censored) rsqrtheta scatter"'

The curve in the qnorm scatter plot is quite linear with some
curvature around the line.

The curvatures around the line could indicate a mixture of models.

Comparing the estimated and the empirical percentiles in the log window suggests
that a good cut-off would be between 5.5 and 6.5. The value 6 is tried.

In the rsqrtheta graph, there is little difference between the lowest and
highest rsquare, so most transformations would do. We choose the value 1.

Run the ros command to assess data with some contamination by outliers.

`". ros yc if yc < 6, censor(censored) scatter"'

Stored results

Scalars

r(r2_a) The adjusted R square from the ROS estimation.

r(rosmean) The estimated mean from ROS.

r(rossd) The estimated standard deviation from ROS.

Matrices

r(ros) Returned cut-off values as specified in option percentiles.

r(rsqrtheta) Returned point values for the rsqrtheta graph.

Author and support

Author:

Niels Henrik Bruun,

Aalborg University Hospital

Support:

niels.henrik.bruun@gmail.com

Development in collaboration with:

Stine Linding Andersen,

Department of Clinical Biochemistry,

Aalborg University Hospital

Nanna Maria Uldall Torp,

Department of Clinical Biochemistry,

Aalborg University Hospital

Peter Astrup Christensen,

Department of Clinical Biochemistry,

Aalborg University Hospital

References

Helsel, Dennis. 2010. "Much Ado About Next to Nothing: Incorporating Nondetects in Science."

The Annals of Occupational Hygiene 54 (April): 257–62. "https://doi.org/10.1093/annhyg/mep092".

Hewett, Paul, and Gary Ganser. 2007. "A Comparison of Several Methods for Analyzing Censored Data."

The Annals of Occupational Hygiene 51 (November): 611–32. "https://doi.org/10.1093/annhyg/mem045".

Hoffmann, Robert G. 1963. "Statistics in the Practice of Medicine." JAMA 185 (11): 864–73.
"https://doi.org/10.1001/jama.1963.03060110068020".

Huston, C., and E. Juarez-Colunga. 2009. "Guidelines for Computing Summary Statistics for Data-Sets Containing Non-Detects."
"https://bvcentre.ca/files/research_reports/08-03GuidanceDocument.pdf".

Jensen, Esther A., Per Hyltoft Petersen, Ole Blaabjerg, Pia Skov Hansen, Thomas H. Brix, and Laszlo Hegedüs. 2006. "Establishment of
Reference Distributions and Decision Values for Thyroid Antibodies Against Thyroid Peroxidase (Tpoab), Thyroglobulin (Tgab) and
the Thyrotropin Receptor (Trab)." Clinical Chemistry and Laboratory Medicine (CCLM) 44 (8): 991–98.
"https://doi.org/doi:10.1515/CCLM.2006.166".