capacity_logreg_main: Estimate channel capacity between discrete input and...

In SLEMI: Statistical Learning Based Estimation of Mutual Information

Estimate channel capacity between discrete input and continuous output

Description

The main wrapping function for basic usage of SLEMI package for estimation of channel capacity. Firstly, data is pre-processed (all arguments are checked, observation with NAs are removed, variables are scaled and centered (if scale=TRUE)). Then basic estimation is carried out and (if testing=TRUE) diagnostic tests are computed. If output directory path is given (output_path is not NULL), graphs visualising the data and the analysis are saved there, together with a compressed output object (as .rds file) with full estimation results.

Usage

capacity_logreg_main(
  dataRaw,
  signal = "input",
  response = NULL,
  output_path = NULL,
  side_variables = NULL,
  formula_string = NULL,
  cc_maxit = 100,
  lr_maxit = 1000,
  MaxNWts = 5000,
  testing = FALSE,
  model_out = TRUE,
  scale = TRUE,
  TestingSeed = 1234,
  testing_cores = 1,
  boot_num = 10,
  boot_prob = 0.8,
  sidevar_num = 10,
  traintest_num = 10,
  partition_trainfrac = 0.6,
  plot_width = 6,
  plot_height = 4,
  data_out = FALSE
)

Arguments

dataRaw	must be a data.frame object
signal	is a character object with names of columns of dataRaw to be treated as channel's input.
response	is a character vector with names of columns of dataRaw to be treated as channel's output
output_path	is the directory in which output will be saved
side_variables	(optional) is a character vector that indicates side variables' columns of data, if NULL no side variables are included
formula_string	(optional) is a character object that includes a formula syntax to use in logistic regression model. If NULL, a standard additive model of response variables is assumed. Only for advanced users.
cc_maxit	is the number of iteration of iterative optimisation of the algorithm to estimate channel capacity. Default is 100.
lr_maxit	is a maximum number of iteration of fitting algorithm of logistic regression. Default is 1000.
MaxNWts	is a maximum acceptable number of weights in logistic regression algorithm. Default is 5000.
testing	is the logical indicating if the testing procedures should be executed
model_out	is the logical indicating if the calculated logistic regression model should be included in output list
scale	is a logical indicating if the response variables should be scaled and centered before fitting logistic regression
TestingSeed	is the seed for random number generator used in testing procedures
testing_cores	- number of cores to be used in parallel computing (via doParallel package)
boot_num	is the number of bootstrap tests to be performed. Default is 10, but it is recommended to use at least 50 for reliable estimates.
boot_prob	is the proportion of initial size of data to be used in bootstrap
sidevar_num	is the number of re-shuffling tests of side variables to be performed. Default is 10, but it is recommended to use at least 50 for reliable estimates.
traintest_num	is the number of overfitting tests to be performed. Default is 10, but it is recommended to use at least 50 for reliable estimates.
partition_trainfrac	is the fraction of data to be used as a training dataset
plot_width	- basic dimensions (width) of plots, in inches
plot_height	- basic dimensions (height) of plots, in inches
data_out	is the logical indicating if the data should be included in output list

Details

In a typical experiment aimed to quantify information flow a given signaling system, input values x_1\leq x_2 \ldots... \leq x_m, ranging from 0 to saturation are considered. Then, for each input level, x_i, n_i observations are collected, which are represented as vectors

y^i_j \sim P(Y|X = x_i)

Within information theory the degree of information transmission is measured as the mutual information

MI(X,Y) = \sum_{i=1}^{m} P(x_i)\int_{R^k} P(y|X = x_i)log_2\frac{P(y|X = x_i)}{P(y)}dy,

where P(y) is the marginal distribution of the output. MI is expressed in bits and 2^{MI} can be interpreted as the number of inputs that the system can resolve on average.

The maximization of mutual information with respect to the input distribution, P(X), defines the information capacity, C. Formally,

C^* = max_{P(X)} MI(X,Y)

Information capacity is expressed in bits and 2^{C^*} can be interpreted as the maximal number of inputs that the system can effectively resolve.

In contrast to existing approaches, instead of estimating, possibly highly dimensional, conditional output distributions P(Y|X =x_i), we propose to estimate the discrete, conditional input distribution, P(x_i |Y = y), which is known to be a simpler problem. Estimation of the MI using estimates of P(x_i |Y = y), denoted here as \hat{P}(x_i|Y = y), is possible as the MI, can be alternatively written as

MI(X,Y) = \sum_{i=1}^{m} P(x_i)\int_{R^k} P(y|X = x_i)log_2\frac{P(x_i|Y = y)}{P(x_i)}dy

The expected value (as in above expression) with respect to distribution P(Y|X = x_i) can be approximated by the average with respect to data

MI(X,Y) \approx \sum_{i=1}^{m} P(x_i)\frac{1}{n_i} \sum_{j=1}^{n_i} P(y|X = x_i)log_2\frac{\hat{P}(x_i|Y = y^i_j)}{P(x_i)}dy

Here, we propose to use logistic regression as \hat{P}(x_i|Y = y^i_j). Specifically,

log\frac{P(x_i |Y = y)}{P(x_m|Y = y)} \approx \alpha_i +\beta_iy

Following this approach, channel capacity can be calculated by optimising MI with respect to the input distribution, P(X). However, this, potentially difficult problem, can be divided into two simpler maximization problems, for which explicit solutions exist. Therefore, channel capacity can be obtained from the two explicit solutions in an iterative procedure known as alternate maximization (similarly as in Blahut-Arimoto algorithm) [1].

Additional parameters: lr_maxit and maxNWts are the same as in definition of multinom function from nnet package. An alternative model formula (using formula_string arguments) should be provided if data are not suitable for description by logistic regression (recommended only for advanced users). Preliminary scaling of data (argument scale) should be used similarly as in other data-driven approaches, e.g. if response variables are comparable, scaling (scale=FALSE) can be omitted, while if they represent different phenomenon (varying by units and/or magnitude) scaling is recommended.

Value

a list with several elements:

• output$regression - confusion matrix of logistic regression predictions

• output$cc - channel capacity in bits

• output$p_opt - optimal probability distribution

• output$model - nnet object describing logistic regression model (if model_out=TRUE)

• output$params - parameters used in algorithm

• output$time - computation time of calculations

• output$testing - a 2- or 4-element output list of testing procedures (if testing=TRUE)

• output$testing_pv - one-sided p-values of testing procedures (if testing=TRUE)

• output$data - raw data used in analysis

References

[1] Csiszar I, Tusnady G, Information geometry and alternating minimization procedures, Statistics & Decisions 1 Supplement 1 (1984), 205–237.

[2] Jetka T, Nienaltowski K, Winarski T, Blonski S, Komorowski M, Information-theoretic analysis of multivariate single-cell signaling responses using SLEMI, PLoS Comput Biol, 15(7): e1007132, 2019, https://doi.org/10.1371/journal.pcbi.1007132.

Examples

tempdata=data_example1
outputCLR1=capacity_logreg_main(dataRaw=tempdata,
signal="signal", response="response",cc_maxit = 10,
formula_string = "signal~response")

tempdata=data_example2
outputCLR2=capacity_logreg_main(dataRaw=tempdata,
signal="signal", response=c("X1","X2"),cc_maxit = 10,
formula_string = "signal~X1+X2") 

#For further details see vignette

SLEMI documentation built on Nov. 20, 2023, 1:06 a.m.