Using cointf from R¶

The R package in r/ is a thin interface to the Python package through reticulate. Every R function calls the Python function of the same name, so results are identical to Python, down to the simulated data for a given seed.

Installation¶

From the root of a clone of the repository:

Rscript -e 'install.packages("reticulate")'   # version 1.41 or later
R CMD INSTALL r

How R finds Python¶

The R package needs a Python environment that has cointf installed. There are two ways to get one.

Automatic (default).
- When the package is loaded, it declares the requirement cointf @ git+https://github.com/aabiddanda/cointf with reticulate::py_require().
- The first time Python is needed, reticulate builds a private environment and installs cointf and its dependencies into it. This needs internet access once and takes about 20 seconds. Later sessions reuse the cached environment.
Your own environment.
- If you installed the Python package with pip install ., point reticulate at that Python before loading the R package.
- The GitHub requirement is then not declared, so local changes to the Python code are what R uses.

Sys.setenv(RETICULATE_PYTHON = "/path/to/venv/bin/python")
library(cointf)

Two options adjust this:

options(cointf.python = "<requirement>"), set before library(cointf), changes the declared requirement, e.g. to another branch (cointf @ git+https://github.com/aabiddanda/cointf@dev) or a local path.
options(cointf.python = NULL) declares nothing. Use it when you select the environment another way, such as reticulate::use_virtualenv().

Loading the R package does not start Python; that happens on the first call. To see which Python is in use, run reticulate::py_config().

Functions¶

R	Python	Purpose
`loglik_fast(data, L, nu, p, kernel, per_meiosis)`	`loglik_fast`	Exact log-likelihood of single meioses
`loglik_quartet(data, L, nu, p, kernel, per_meiosis)`	`loglik_quartet`	Log-likelihood of phase-unknown (quartet) data
`meioses(data, L)`	`Meioses`	Prepare data once for repeated evaluation
`chebyshev_kernel(...)`	`ChebyshevKernel`	Faster approximate kernel, cached per `nu`
`simulate_meioses(L, nu, p, n, seed)`	`simulate_meioses`	Simulate meioses
`simulate_quartets(L, nu, p, n, seed)`	`simulate_quartets`	Simulate quartet data
`meiosis_list(position, meiosis)`	(R only)	Convert long-format data to the list format
`cointf_python()`	the `cointf` module	Direct access to anything without a wrapper

Each function has a help page, e.g. ?loglik_fast.

Data¶

Crossover data are a list with one numeric vector of positions (in Morgans) per meiosis. Use numeric(0) for a meiosis without crossovers.

data <- list(c(0.31, 1.12), 0.70, numeric(0))
loglik_fast(data, 1.5, nu = 4, p = 0.05)
loglik_fast(data, 1.5, nu = 4, p = 0.05, per_meiosis = TRUE)

Data often come as one row per crossover. meiosis_list() converts that format. A meiosis without crossovers is kept by giving it one row with NA position:

df <- data.frame(id = c("a", "a", "b", "c"), pos = c(0.2, 0.9, NA, 0.5))
data <- meiosis_list(df$pos, df$id)   # list(a = c(0.2, 0.9), b = numeric(0), c = 0.5)

L, nu and p each accept one value or one value per meiosis, as in Python.

Repeated evaluation¶

Each call from an R list sends the data to Python and builds the internal data structure (Meioses). For 20,000 meioses that takes about 80 ms, most of it in building the structure. When the same data are evaluated many times, as in an optimiser, build the structure once with meioses() and pass it with L = NULL. Each evaluation then costs only the likelihood itself (a few ms for 20,000 meioses with the Chebyshev kernel).

mz <- meioses(data, 1.5)
kernel <- chebyshev_kernel()          # keep one kernel so its cache is reused
loglik_fast(mz, NULL, nu = 4, p = 0.05, kernel = kernel)
loglik_quartet(mz, NULL, nu = 4, p = 0.05, kernel = kernel)   # same object works for quartets

The object returned by meioses() and the kernel are Python objects. They stay in Python between calls; only numbers cross back into R.

Example: maximum-likelihood fit¶

The same fit as the Python example, working on \((\log\nu, \operatorname{logit} p)\):

library(cointf)

data <- simulate_meioses(L = 1.5, nu = 4, p = 0.05, n = 2000, seed = 1)
mz <- meioses(data, 1.5)
kernel <- chebyshev_kernel()

negloglik <- function(theta) {
  -loglik_fast(mz, NULL, nu = exp(theta[1]), p = plogis(theta[2]), kernel = kernel)
}

fit <- optim(c(log(2), -2), negloglik, hessian = TRUE)
se <- sqrt(diag(solve(fit$hessian)))
est <- fit$par
ci <- rbind(est - 1.96 * se, est + 1.96 * se)

data.frame(
  estimate = c(nu = exp(est[1]), p = plogis(est[2])),
  lower = c(exp(ci[1, 1]), plogis(ci[1, 2])),
  upper = c(exp(ci[2, 1]), plogis(ci[2, 2]))
)
#>      estimate      lower     upper
#> nu 3.77266971 3.30427871 4.3074565
#> p  0.04922224 0.02937439 0.0813569

The fit takes about 0.06 s and agrees with the Python notebook to every printed digit (the data are identical for the same seed). For phase-unknown data, replace loglik_fast with loglik_quartet.

Testing¶

export RETICULATE_PYTHON=$(which python)                      # Python with your local install
Rscript -e 'testthat::test_local("r")'                        # full R suite
COINTF_TEST_INSTALL=true Rscript -e 'testthat::test_local("r")'   # also test the GitHub install

The R tests cover:

the Python reference values, the quartet closed form and brute force, and the Chebyshev kernel;
that results are identical to calling Python directly;
conversion of awkward inputs (empty, named, integer, unsorted, NA/Inf) and per-meiosis parameters;
Python objects surviving round trips through R;
a bound on conversion time;
in fresh R sessions, that loading the package does not start Python and that RETICULATE_PYTHON is honoured without a warning.