Implementation Overview

1 Introduction

This vignette discusses the implementation of the package DistrFree. At its core is the constructor function DistrFree() which currently represents the distribution of a model of the form

\[ F_Y(y \given G = k) = F_Z(F_Z^{-1}(F_Y(y \given G = 1)) - \delta_k). \tag{1}\]

We have $K$ treatment groups $G \in \{ 1, \dots, K \}$, $K \geq 2$ for an at least ordered outcome $Y \in \mathcal{Y} = \{v_1, v_2, \dots, v_C\}$, where $C$ is the number of unique support points. We refer to the first treatment $G = 1$ as control with $\delta_1 := 0$. The model is parameterized through some absolutely continuous CDF $F_Z$ with log-concave density $f_Z = F'_Z$ and corresponding derivative $f'_Z$. The link function $F_Z^{-1}$ defines the interpretation of the shift parameters $\delta_k$ (see Table 1), and positive values of the parameters $\delta_k$ correspond to stochastically larger response values in group $k$ when compared to the control group (Hothorn and Hornik 2026).

Table 1: Interpretation of the shift parameters $\delta_k$ depending on the choice of $F_Z$, i.e. the link argument in the constructor function.

$F_Z$	`link`	Interpretation of $\delta_k$
$\exp(-\exp(-z))$	`"loglog"`	log-reverse time hazard ratio
$1-\exp(-\exp(z))$	`"cloglog"`	log-hazard ratio
$\expit(z)$	`"logit"`	log-odds ratio
$\Phi(z)$	`"probit"`	generalized Cohen’s $d$

The constructor function DistrFree() takes three arguments in order to specify distributions from the model class in Equation 1:

baseline represents the CDF of the control group, i.e. $F_Y(y \given G = 1)$. If the user specifies baseline through a numeric vector, a discrete uniform distribution is created, see Listing 1 as an example. The user can also provide some data.frame with two columns, of which one is either called weights or probability, in order to have more control over the baseline distribution.
location is a numeric vector defining the shift effects $\delta_k$ on the latent scale defined by $F_Z$. Each element defines one distribution which share their baseline and link. The DistrFree object created in Listing 1 contains the control group ($\delta_1 = 0$) and a second group with $\delta_2 = 2$.
link represents the link function $F_Z^{-1}$ and defaults to "logit".

Listing 1: Construction of a DistrFree object.

d <- DistrFree(location = c(0, 2), baseline = seq_len(100), link = "logit")
d

[1] "DistrFree(location = 0)" "DistrFree(location = 2)"

2 Methods

2.1 Compute the Mean, Variance, Skewness and Kurtosis

mean(d)

[1] 50.50 79.94

variance(d)

[1] 833.2 499.3

skewness(d)

[1]  3.582e-18 -1.536e+00

kurtosis(d)

[1] -1.200  1.741

2.2 Compute the Cumulative Distribution Function

For some arbitrary $y \in \mathbb{R}$, the cumulative distribution function (CDF) of $Y \in \mathcal{Y}$ for treatment group $G = k$ is given through \[ \begin{aligned} F_{Y}(y \given G = k) &= \mathbb{P}_Y(Y \in \{ v \in \mathcal{Y} : v \leq y \}) \\ &= F_Z(F_Z^{-1}(F_{Y}(y \given G = 1)) - \delta_k) \end{aligned}, \]

as denoted in Equation 1. Consequently, the CDF is a right-continuous step function jumping at the support points $\{v_1, v_2, \dots, v_C\}$ and constant in between, see Listing 2 and Figure 1 below.

Listing 2: Compute the Cumulative Distribution Function for a DistrFree object.

cdf(d, c(3, 7, 12, 12.5))

          p_3     p_7    p_12  p_12.5
[1,] 0.030000 0.07000 0.12000 0.12000
[2,] 0.004168 0.01008 0.01812 0.01812

Figure 1: Cumulative distribution functions of the control and treatment group.

2.3 Compute the Probability Mass Function

Since $Y \in \mathcal{Y} = \{v_1, v_2, \dots, v_C\}$ with $v_1 < v_2 < \dots < v_C$ is discrete, values not in $\mathcal{Y}$ have zero density. The function pdf(<DistrFree>, x, ...) consequently computes the probability mass for a given value x¹.

The density for such a discrete response is given through \[ f_Y(v_c \given G = k) = F_Y(v_c \given G = k) - F_Y(v_{c-1} \given G = k), \quad c = 1, \dots, C, \]

where $F_Y(v_0) := 0$ (Hothorn et al. 2018). As such, we can compute the density at any existing support point, whereas a value not in the sample space returns zero.

pdf(d, c(1, 25, 50.5))

          d_1     d_25 d_50.5
[1,] 0.010000 0.010000      0
[2,] 0.001365 0.002179      0

We can plot both densities of the DistrFree object generated in Listing 1. Figure 2 illustrates how the probability masses of $\delta_2 = 2$ are shifted towards larger values compared to the control group ($\delta_1 = 0$).

Figure 2: Probability mass functions of the control and treatment group.

2.4 Compute the Quantile Function

The function quantile(<DistrFree>, probs, ...) returns the quantile $Q_Y(p \given G = k) = F_Y^{-1}(p \given G = k) = \inf\{ y \in \mathcal{Y} : F_Y(y \given G = k) \geq p \}$ for $p \in [0, 1]$.

quantile(d, probs = c(0, 0.25, 0.75, 1.1, NA))

Warning in .qfun(p = at, d = d, ...): NaNs produced

     q_0 q_0.25 q_0.75 q_1.1 q_NA
[1,]   1     25     75   NaN   NA
[2,]   1     72     96   NaN   NA

We can furthermore plot the quantile functions $Q_Y(p \given G = k)$ across $p \in [0, 1]$ in Figure 3.

Figure 3: Quantile functions of the control and treatment group.

2.4.1 Rounding Error

As in the function q_DISCRETE_BODY() defined in the R source file src/nmath/qDiscrete_search.h (R Core Team 2026), quantile.DistrFree() contains fuzz to ensure left-continuity at the price of possible errors in the upper tail.

The issue with rounding error is visible in the following simple example. For baseline = seq_len(3), the sample space is $\mathcal{Y} = \{v_1, v_2, v_3\} = \{1, 2, 3\}$. We compute $F_Y(y \given G = 2)$ as defined in Equation 1, where $F_Z = \expit$, $F_Z^{-1} = \logit$ and $\delta_2 = \log(2)$.

d0 <- DistrFree(location = log(2), baseline = seq_len(3), link = "logit")

link <- attr(d0, "link")
support <- attr(d0, "support")
baseline <- attr(d0, "baseline")

h <- link$linkfun(baseline(support))
cdf <- link$linkinv(h - d0$location)

By choosing $\delta_2 = \log(2)$, the odds of being in category $v_1$ are halved when moving from the control group to $G = 2$, i.e.

\[ \begin{aligned} F_Y(v_1 \given G = 2) &= \expit(\logit(1 / 3) - \log(2)) \\ &= \expit(\log(1 / 4)) \\ &= 1 / 5 \end{aligned} \]

As such we would expect that $F_Y^{-1}(0.2 \given G = 2) = v_1 = 1$, but the output disagrees.

prb <- 0.2
support[which(cdf >= prb)[1L]]

[1] 2

Reason for this disagreement is some rounding error which is present in the computed CDF.

print(cdf, digits = 20)

[1] 0.19999999999999995559 0.50000000000000000000 1.00000000000000000000

If we add some fuzz, currently hardcoded to 64 * .Machine$double.eps, the rounding error is cancelled in this scenario, and we arrive at the expected $F_Y^{-1}(0.2 \given G = 2) = v_1 = 1$.

eps <- 64 * .Machine$double.eps
support[which(cdf >= (prb * (1 - eps)))[1L]]

[1] 1

2.5 Draw Random Samples

As stated in Hothorn and Hornik (2026), we know with Equation 1 that for an absolutely continuous random variable $Y$ \[ U = F_Y(Y \given G = k) = F_Z(F_Z^{-1}(F_Y(Y \given G = 1)) - \delta_k), \quad k = 2, \dots, K \]

follows a standard uniform distribution on the unit interval, allowing us to sample from the distribution of $Y$ using \[ F_Y^{-1}(F_Z(F_Z^{-1}(U) + \delta_k) \given G = 1). \]

Assuming a uniform distribution for $F_Y$ in the control group, it is enough to draw samples from $F_Z(F_Z^{-1}(U) + \delta_k)$. Discretization follows by transforming the drawn samples of all groups by the same quantile function $F_Y^{-1}(p \given G = 1) = \inf\{ y \in \mathcal{Y} : F_Y(y \given G = 1) \geq p \}$ post-hoc, which does not change the shift effects.

The implementation of random(<DistrFree>, n, ...) follows the logic of the internal function stats:::.rfree1way() closely (R Core Team 2026).

set.seed(29)
random(d, n = 10)

     r_1 r_2 r_3 r_4 r_5 r_6 r_7 r_8 r_9 r_10
[1,]  10  11  59  83  12  99  31  18  67   87
[2,]  71  79  44  98  70  83  93  98  81   81

References

Hothorn, Torsten, and Kurt Hornik. 2026. Free1way.docreg: Additional Documentation and Regression Tests for ’Stats::free1way()’. https://codeberg.org/thothorn/coin/.

Hothorn, Torsten, Lisa Möst, and Peter Bühlmann. 2018. “Most Likely Transformations.” Scandinavian Journal of Statistics 45 (1): 110–34. https://doi.org/10.1111/sjos.12291.

R Core Team. 2026. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. https://doi.org/10.32614/R.manuals.

Footnotes

Recall that the function pmf() exists, though, it is simply an alias which internally calls pdf().↩︎

\(F_Z\)	`link`	Interpretation of \(\delta_k\)
\(\exp(-\exp(-z))\)	`"loglog"`	log-reverse time hazard ratio
\(1-\exp(-\exp(z))\)	`"cloglog"`	log-hazard ratio
\(\expit(z)\)	`"logit"`	log-odds ratio
\(\Phi(z)\)	`"probit"`	generalized Cohen’s \(d\)