README.Rmd

---
output: 
  github_document:
    html_preview: false
---

```{r, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-"
)
```

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, logo, include=FALSE}

library(hexSticker)
library(ggplot2)

dbinorm <- function(x, p, dif) {
  p * dnorm(x, 0, 1) + (1 - p) * dnorm(x, dif, 1)
}
dbinorm_scaled <- function(x, scaling, ...) {
  dbinorm(x, ...) * scaling
}
dnorm_scaled <- function(x, scaling, ...) {
  dnorm(x, ...) * scaling
}

dratio <- function(x, p, dif, mu, sd) {
  dbinorm(x, p, dif) / dnorm(x, mu, sd)
}

hexSticker::sticker(
    ggplot() +
        stat_function(aes(col = "frac(italic(P[n](x)), italic(P[d](x)))"),
                      fun = dratio, args = list(p = 0.4, dif = 3, mu = 3, sd = 2),
                      alpha = 0.9, size = 0.9) +
        stat_function(fun = dbinorm_scaled, args = list(p=0.4, dif = 3, scaling = 10),
                      alpha = 0.3, col = "white", size = 0.3, linetype = "dashed") +
        stat_function(fun = dnorm_scaled, args = list(scaling = 10, mean = 3, sd = 2),
                      alpha = 0.3, col = "white", size = 0.3, linetype = "dashed") +

        theme_void() +
        hexSticker::theme_transparent() +
        scale_color_manual(name = NULL, values = c("#d8e5e8"), labels = scales::parse_format()) +
        theme(panel.background = element_rect(fill = 'transparent'),
              plot.background = element_rect(fill = "transparent"),
              legend.position = c(1.02, 0.53),
              legend.text = element_text(colour = "white", size = 24)) +
        guides(color = guide_legend(override.aes = list(linetype = 0))) +
        xlim(-5, 8.6),
    package = "densityratio",
    p_size = 22, s_x = 1, s_y = 1, s_width = 1.4, s_height = 1.4,
    h_fill = "#034e57", h_color = "#93effa", p_color = "#de0277",
    p_y = 1.4,
    filename = "man/figures/logo.png"
)
```

# densityratio <img src="man/figures/logo.png" align="right" alt="densratio logo" width="130"/>

<!-- badges: start -->
[![R-CMD-check](https://github.com/thomvolker/densityratio/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/thomvolker/densityratio/actions/workflows/R-CMD-check.yaml)
[![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
[![densityratio status badge](https://thomvolker.r-universe.dev/badges/densityratio)](https://thomvolker.r-universe.dev/densityratio)
[![DOI](https://zenodo.org/badge/613441108.svg)](https://zenodo.org/badge/latestdoi/613441108)

<!-- badges: end -->

## Overview

This package provides functionality to directly estimate a density ratio
$$r(x) = \frac{p_\text{nu}(x)}{p_{\text{de}}(x)},$$
without estimating the numerator and denominator density separately. Density ratio estimation serves many purposes, for example, prediction, outlier detection, change-point detection in time-series, importance weighting under domain adaptation (i.e., sample selection bias) and evaluation of synthetic data utility. The key idea is that differences between data distributions can be captured in their density ratio, which is estimated over the entire multivariate space of the data. Subsequently, the density ratio values can be used to summarize the dissimilarity between the two distributions in a discrepancy measure.

```{r, densities, echo=FALSE,message=FALSE,warning=FALSE,fig.height=6,fig.width=10,out.width="15cm",fig.dpi=500,fig.format="png",fig.align="center"}

set.seed(1)
library(ggplot2)
library(patchwork)

dbinorm <- function(x, p, dif) {
  p * dnorm(x, 0, 1) + (1 - p) * dnorm(x, dif, 1)
}

dratio <- function(x, p, dif, mu, sd) {
  dbinorm(x, p, dif) / dnorm(x, mu, sd)
}

fit <- densityratio::ulsif(
  df_numerator = densityratio::numerator_data$x5, 
  df_denominator = densityratio::denominator_data$x5,
  centers = c(densityratio::numerator_data$x5, densityratio::denominator_data$x5),
  parallel = TRUE,
  nthreads = 18
)

ggplot() +
  stat_function(aes(col = factor(1, labels = "italic(P[n](x))")), 
                fun = dbinorm, args = list(p = 0.4, dif = 3),
                size = 1) +
  stat_function(aes(col = factor(2, labels = "italic(P[d](x))")), 
                fun = dnorm, args = list(mean = 3, sd = 2),
                size = 1) +
  xlim(-5, 8.6) +
  scale_color_manual(name = NULL, 
                     values = c("#93effa", "#034e57"), 
                     labels = scales::parse_format()) +
  theme_classic() +
  hexSticker::theme_transparent() +
  theme(legend.position = c(0.9, 0.9),
        text = element_text(size = 20),
        panel.background = element_rect(fill = 'transparent'),
        plot.background = element_rect(fill = "transparent")) +
  xlab(expression(italic(x))) +
  ylab(NULL) +
ggplot() +
  stat_function(aes(col = "italic(r(x))"), 
                fun = dratio, args = list(p = 0.4, dif = 3, mu = 3, sd = 2),
                size = 1) +
  geom_line(aes(x = seq(-5, 8.6, length.out = 100),
                y = predict(fit, seq(-5, 8.6, length.out = 100)),
                col = "italic(hat(r)(x))"),
            size = 1) +
  xlim(-5, 8.6) +
  theme_classic() +
  hexSticker::theme_transparent() +
  scale_color_manual(name = NULL, 
                     values = c("pink", "#de0277"),
                     labels = scales::parse_format()) +
  theme(legend.position = c(0.9, 0.9),
        text = element_text(size = 20),
        panel.background = element_rect(fill = 'transparent'),
        plot.background = element_rect(fill = "transparent")) +
  xlab(expression(italic(x))) +
  ylab(NULL)
```


### Features

* __Fast__: Computationally intensive code is executed in `C++` using `Rcpp` and `RcppArmadillo`.
* __Automatic__: Good default hyperparameters that can be optimized in cross-validation (we do recommend understanding those parameters before using `densityratio` in practice).
* __Complete__: Several density ratio estimation methods, such as unconstrained least-squares importance fitting (`ulsif()`), Kullback-Leibler importance estimation procedure (`kliep()`), ratio of estimated densities (`naive()`), ratio of estimated densities after dimension reduction (`naivesubspace()`), and least-squares heterodistributional subspace search (`lhss()`; experimental).
* __User-friendly__: Simple user interface, default `predict()`, `print()` and `summary()` functions for all density ratio estimation methods; built-in data sets for quick testing.


## Installation

You can install the development version of `densityratio` from [R-universe](https://r-universe.dev/search) with:

``` r
install.packages('densityratio', repos = 'https://thomvolker.r-universe.dev')
```

## Usage

The package contains several functions to estimate the density ratio between the numerator data and the denominator data. To illustrate the functionality, we make use of the in-built simulated data sets `numerator_data` and `denominator_data`, that both consist of the same five variables.

### Minimal example

```{r, example-univ}

library(densityratio)

head(numerator_data)

fit  <- ulsif(
  df_numerator = numerator_data$x5, 
  df_denominator = denominator_data$x5, 
  nsigma = 5, 
  nlambda = 5
)

class(fit)
```

We can ask for the `summary()` of the estimated density ratio object, that contains the optimal kernel weights (optimized using cross-validation) and a measure of discrepancy between the numerator and denominator densities.

```{r, summary-univ}
summary(fit)
```

To formally evaluate whether the numerator and denominator densities differ significantly, you can perform a two-sample homogeneity test as follows.

```{r, summary-univ-test}
summary(fit, test = TRUE)
```

The probability that numerator and denominator samples share a common data generating mechanism is very small.

The `ulsif`-object also contains the (hyper-)parameters used in estimating the density ratio, such as the centers used in constructing the Gaussian kernels (`fit$centers`), the different bandwidth parameters (`fit$sigma`) and the regularization parameters (`fit$lambda`). 
Using these variables, we can obtain the estimated density ratio using `predict()`.


```{r, plot-univ, fig.height=6, fig.width=10, out.width="15cm", fig.dpi=500, fig.format="png", fig.align="center"}

# obtain predictions for the numerator samples

newx5 <- seq(from = -3, to = 6, by = 0.05)
pred  <- predict(fit, newdata = newx5)

ggplot() +
  geom_point(aes(x = newx5, y = pred, col = "ulsif estimates")) +
  stat_function(mapping = aes(col = "True density ratio"), 
                fun = dratio, 
                args = list(p = 0.4, dif = 3, mu = 3, sd = 2),
                linewidth = 1) +
  theme_classic() +
  scale_color_manual(name = NULL, values = c("#de0277", "purple")) +
  theme(legend.position.inside = c(0.8, 0.9),
        text = element_text(size = 20))
```

### Scaling

By default, all functions in the `densityratio` package standardize the data to the numerator means and standard deviations. This is done to ensure that the importance of each variable in the kernel estimates is not dependent on the scale of the data. By setting `scale = "denominator"` one can scale the data to the means and standard deviations of the denominator data, and by setting `scale = FALSE` the data remains on the original scale.


### Categorical data

All of the functions in the `densityratio` package accept categorical variables types. However, note that internally, these variables are one-hot encoded, which can lead to a high-dimensional feature-space. 


```{r, cat-data}

summary(numerator_data$x1)
summary(denominator_data$x1)

fit_cat <- ulsif(
  df_numerator = numerator_data$x1, 
  df_denominator = denominator_data$x1
)

aggregate(
  predict(fit_cat) ~ numerator_data$x1,
  FUN = unique
)

table(numerator_data$x1) / table(denominator_data$x1)
```

This transformation can give a reasonable estimate of the ratio of proportions in the different data sets (although there is some regularization applied such that the estimated odds are closer to one than seen in the real data).


### Full data example

After transforming all variables to numeric variables, it is possible to calculate the density ratio over the entire multivariate space of the data.

```{r, full-data-example}
fit_all <- ulsif(
  df_numerator = numerator_data,
  df_denominator = denominator_data
)

summary(fit_all, test = TRUE, parallel = TRUE)
```


### Other density ratio estimation functions

Besides `ulsif()`, the package contains several other functions to estimate a density ratio. 

* `naive()` estimates the numerator and denominator densities separately, and subsequently takes there ratio. 
* `kliep()` estimates the density ratio directly through the Kullback-Leibler importance estimation procedure.
<!-- * `kmm()` estimates the density ratio for the denominator sample points only through kernel mean matching.  -->

```{r, plot-methods, warning=FALSE, fig.height=6, fig.width=10, out.width="15cm", fig.dpi=500, fig.format="png", fig.align="center"}

fit_naive <- naive(
  df_numerator = numerator_data$x5, 
  df_denominator = denominator_data$x5
)

fit_kliep <- kliep(
  df_numerator = numerator_data$x5, 
  df_denominator = denominator_data$x5
)

pred_naive <- predict(fit_naive, newdata = newx5)
pred_kliep <- predict(fit_kliep, newdata = newx5)

ggplot(data = NULL, aes(x = newx5)) +
  geom_point(aes(y = pred, col = "ulsif estimates")) +
  geom_point(aes(y = pred_naive, col = "naive estimates")) +
  geom_point(aes(y = pred_kliep, col = "kliep estimates")) +
  stat_function(aes(x = NULL, col = "True density ratio"), 
                fun = dratio, args = list(p = 0.4, dif = 3, mu = 3, sd = 2),
                linewidth = 1) +
  theme_classic() +
  scale_color_manual(name = NULL, values = c("pink", "#512970","#de0277", "purple")) +
  theme(legend.position.inside = c(0.8, 0.9),
        text = element_text(size = 20))
```

The figure directly shows that `ulsif()` and `kliep()` come rather close to the true density ratio function in this example, and outperform the `naive()` solution. 

## Contributions

This package is still in development, and I'll be happy to take feedback and suggestions. Please submit these through [GitHub Issues](https://github.com/thomvolker/densityratio/issues). 

## Resources


__Books__

* General information about the density ratio estimation framework: Sugiyama, Suzuki and Kanamori (2012). [Density Ratio Estimation in Machine Learning](https://www.cambridge.org/core/books/density-ratio-estimation-in-machine-learning/BCBEA6AEAADD66569B1E85DDDEAA7648)


__Papers__

* Density ratio estimation for the evaluation of the utility of synthetic data: Volker, De Wolf and Van Kesteren (2023). [Assessing the utility of synthetic data: A density ratio perspective](https://unece.org/statistics/documents/2023/08/working-documents/assessing-utility-synthetic-data-density-ratio)

## How to cite

Volker, T.B. (2023). densityratio: Distribution comparison through density ratio estimation. <https://doi.org/10.5281/zenodo.8307819>