--- title: "Implementation Overview" format: html: embed-resources: true html-math-method: mathjax toc: true number-sections: true link-citations: true bibliography: ../inst/REFERENCES.bib vignette: > %\VignetteIndexEntry{Implementation Overview} %\VignetteEngine{quarto::format} %\VignetteDepends{DistrFree,distributions3} %\VignetteKeywords{distribution-free} %\VignettePackage{DistrFree} %\VignetteEncoding{UTF-8} --- ```{r preliminaries, echo = FALSE, message = FALSE} suppressMessages(library("DistrFree")) options(digits = 4) knitr::opts_chunk$set(fig.align = "center") ``` ```{css, echo = FALSE, eval = FALSE} p { text-align: justify; } ``` {{< include _macros.qmd >}} ## 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). $$ {#eq-model-definition} 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 @tbl-choice-linkinv), and positive values of the parameters $\delta_k$ correspond to stochastically larger response values in group $k$ when compared to the control group [@free1way2026]. | $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$ | : Interpretation of the shift parameters $\delta_k$ depending on the choice of $F_Z$, *i.e.* the `link` argument in the constructor function. {#tbl-choice-linkinv} The constructor function `DistrFree()` takes three arguments in order to specify distributions from the model class in @eq-model-definition: - `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 @lst-constructor 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 @lst-constructor 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"`. ```{r constructor, echo = TRUE, message = TRUE} #| lst-label: lst-constructor #| lst-cap: Construction of a `DistrFree` object. d <- DistrFree(location = c(0, 2), baseline = seq_len(100), link = "logit") d ``` ## Methods ### Compute the Mean, Variance, Skewness and Kurtosis ```{r moments, echo = TRUE, message = TRUE} mean(d) variance(d) skewness(d) kurtosis(d) ``` ### 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 @eq-model-definition. 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 @lst-cdf and @fig-cdf below. ```{r cdf, echo = TRUE, message = TRUE} #| lst-label: lst-cdf #| lst-cap: Compute the Cumulative Distribution Function for a `DistrFree` object. cdf(d, c(3, 7, 12, 12.5)) ``` ```{r, echo = FALSE} #| label: fig-cdf #| fig-cap: "Cumulative distribution functions of the control and treatment group." GRID <- seq(-1, 101, 1) CDF <- t(cdf(d, GRID)) graphics::matplot(GRID, CDF, type = "s", xlab = "y", ylab = "F(y)", xlim = c(0, 100), lty = 1, main = "Cumulative Distribution Function") legend("topleft", legend = paste("delta =", d$location), col = seq_along(d$location), bty = "n", lty = 1) ``` ### 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(, x, ...)` consequently computes the probability mass for a given value `x`^[Recall that the function `pmf()` exists, though, it is simply an alias which internally calls `pdf()`.]. 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$ [@Hothorn2018]. As such, we can compute the density at any existing support point, whereas a value not in the sample space returns zero. ```{r pdf, echo = TRUE, message = TRUE} pdf(d, c(1, 25, 50.5)) ``` We can plot both densities of the `DistrFree` object generated in @lst-constructor. @fig-pmf illustrates how the probability masses of $\delta_2 = `r d$location[2]`$ are shifted towards larger values compared to the control group ($\delta_1 = `r d$location[1]`$). ```{r, echo = FALSE} #| label: fig-pmf #| fig-cap: "Probability mass functions of the control and treatment group." GRID <- seq_len(100) PMF <- t(pdf(d, GRID)) graphics::matplot(GRID, PMF, type = "h", xlab = "y", ylab = "f(y)", lty = 1, main = "Probability Mass Function") graphics::matpoints(GRID, PMF, pch = 20, cex = 0.5) legend("topleft", legend = paste("delta =", d$location), col = seq_along(d$location), lty = 1, pch = 20, bty = "n") ``` ### Compute the Quantile Function The function `quantile(, 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]$. ```{r quantile, echo = TRUE, message = TRUE} quantile(d, probs = c(0, 0.25, 0.75, 1.1, NA)) ``` We can furthermore plot the quantile functions $Q_Y(p \given G = k)$ across $p \in [0, 1]$ in @fig-quantile. ```{r, echo = FALSE} #| label: fig-quantile #| fig-cap: "Quantile functions of the control and treatment group." GRID <- seq(0, 1, 0.01) VAL <- t(quantile(d, GRID)) graphics::matplot(GRID, VAL, type = "S", xlab = "p", ylab = "Q(p)", lty = 1, main = "Quantile Function") legend("topleft", legend = paste("delta =", d$location), col = seq_along(d$location), lty = 1, bty = "n") ``` #### Rounding Error As in the function `q_DISCRETE_BODY()` defined in the R source file `src/nmath/qDiscrete_search.h` [@RCoreTeam], `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 @eq-model-definition, where $F_Z = \expit$, $F_Z^{-1} = \logit$ and $\delta_2 = \log(2)$. ```{r quantile_error_01, echo = TRUE, message = TRUE} 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. ```{r quantile_error_02, echo = TRUE, message = TRUE} prb <- 0.2 support[which(cdf >= prb)[1L]] ``` Reason for this disagreement is some rounding error which is present in the computed CDF. ```{r quantile_error_03, echo = TRUE, message = TRUE} print(cdf, digits = 20) ``` 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$. ```{r quantile_error_04, echo = TRUE, message = TRUE} eps <- 64 * .Machine$double.eps support[which(cdf >= (prb * (1 - eps)))[1L]] ``` ### Draw Random Samples As stated in @free1way2026, we know with @eq-model-definition 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(, n, ...)` follows the logic of the internal function `stats:::.rfree1way()` closely [@RCoreTeam]. ```{r random, echo = TRUE, message = TRUE} set.seed(29) random(d, n = 10) ``` ## References {.unnumbered} ::: {#refs} :::