Statistical Learning based Estimation of the Mutual Information (SLEMI) - R package

User Manual

T. Jetka1, K. Nienałtowski, M. Komorowski

22 February 2021

Abstract

The package SLEMI is designed to estimate channel capacity between finite state input and multidimensional continuous output from experimental data. For efficient computations, it uses an iterative algorithm based on logistic regression. In addition, functions to estimate mutual information and calculate probabilities of correct discrimination between a pair of input values are implemented. The method is published in PLOS Computational Biology (Jetka et al. 2019).

1 Preliminaries

1.1 Requirements - Hardware

1.2 Requirements - Software

The main software requirement is the installation of the R environment (version: >= 3.2), which can be downloaded from R project website and is distributed for all common operating systems. We tested the package in R environment installed on Windows 7, 10; Mac OS X 10.11 - 10.13 and Ubuntu 18.04 with no significant differences in the performance. The use of a dedicated Integrated development environment (IDE), e.g. RStudio is recommended.

Apart from a base installation of R, SLEMI requires the following R packages:

  1. for installation
  1. for estimation
  1. for visualisation
  1. for data handling

Each of the above packages can be installed by executing

install.packages("name_of_a_package")

in the R console.

Importantly, during installation availability of the above packages will be verified and missing packages will be automatically installed.

1.3 Installation

The package can be directly installed from GitHub. For installation, open RStudio (or base R) and run following commands in the R console

install.packages("devtools") # run if 'devtools' is not installed
library(devtools)
install_github("sysbiosig/SLEMI")

Are required packages not found, they will be installed automatically.

1.4 Citing and support

The package implements methods published in PLOS Computational Biology, please cite:

Jetka T, Nienałtowski K, Winarski T, Błoński S, Komorowski M (2019) Information-theoretic analysis of multivariate single-cell signaling responses. PLoS Comput Biol 15(7): e1007132. https://doi.org/10.1371/journal.pcbi.1007132

All problems, issues and bugs can be reported here:

https://github.com/sysbiosig/SLEMI/issues

or directly via e-mail: t.jetka a t gmail.com.

2 Structure of the package

The three functions listed below constitute the key wrapper (interface) functions of the package.

  1. mi_logreg_main() enables calculation of the mutual information

  2. capacity_logreg_main() enables calculation of the information capacity

  3. prob_discr_pairwise() serves to calculate probabilities of correct discrimination between pairs of input values

The function capacity_logreg_main() triggers

  1. preprocessing of the data

  2. estimation of channel capacity

  3. running diagnostic procedures

  4. visualisation.

Each of the above steps is implemented within auxiliary functions as presented in the Figure 1 below.

The algorithm to compute the information capacity is implemented within the function capacity_logreg_algorithm(), which uses logistic regression from the nnet package.

Diagnostic procedures (significance and uncertainties of estimates) are provided in an internal function capacity_logreg_testing(). These are based on data bootstrapping and overfitting test.

For visualization, a set of graphs is created by an internal function capacity_output_graphs() and saved in a specified directory. In addition, capacity_logreg_main() returns a list with capacity estimates, optimal input probability distribution, diagnostic measures and other summary information about the analysis.

The function mi_logreg_main() serves to calculate the mutual information. It initiates similar steps as the function capacity_logreg_main() but without performing the optimization of the distribution of the input. Instead, it requires the input distribution to be specified by the user as a function’s argument.

Logistic regression and Monte Carlo methods, following an analogous algorithm as within the capacity_logreg_algorithm() function, are combined to estimate mutual information within a function mi_logreg_algorithm(). Visualisation and diagnostics are carried out by the same set of auxiliary functions as for channel capacity (internal functions capacity_output_graphs() and capacity_logreg_testing()).

The prob_discr_pairwise() allows to estimate probabilities of correct discrimination between two different values of the input. It implements estimation of probabilities of correct classification by logistic regression (from nnet package) for each pair of input values. The probabilities of correct discrimination are visualized with a graph composed of pie charts.

Main function to estimate channel capacity

2.1 Formulation of the problem

SLEMI package is designed to estimate information-theoretic measures between a discrete-valued input, \(X\), and multivariate, continuous output, \(Y\). 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. For details regarding information theory or its application in systems biology please see Methods section and Supplementary Information of the corresponding paper (Jetka et al. 2019).

2.2 Input data

Functions mi_logreg_main(), capacity_logreg_main(), prob_discr_pairwise() require data in the form of the object data.frame with a specific structure of rows and columns. Responses \(y^i_j\) are assumed to be measured for a finite set of stimuli levels \(x_1,x_2,\ldots,x_m\). The responses \(y^i_j\) can be multidimensional. Usually, experimental dataset is represented as a table with rows and columns organized as shown in Figure 2.

Therefore, the input data frame is expected to have the form represented by the above table, which can be formally described by the following conditions

Conceptual representation of a generic experimental dataset needed for quantifying information transmission of a channel

An example of the input data.frame, which contains the measurements of the NfkB system presented in the MP is available within the package under the variable data_nfkb. It has the following format

signal response_0 response_3 response_21
1 0ng 0.3840744 0.4252835 0.5643459
2 0ng 0.4709216 0.5777821 0.2962640
3 0ng 0.4274474 0.6696011 0.5696369
10001 8ng 0.3120216 0.3475484 9.7036535
10002 8ng 0.2544961 0.6611051 8.1628482
10003 8ng 0.1807391 0.4336810 5.3928484
11540 100ng 1.3534083 3.0158004 6.8983046
11541 100ng 1.7007936 2.2224497 2.8677178
11542 100ng 0.1997087 0.2886905 5.8193494

where each row represents measurements of a single-cell, the column named signal specifies the level of stimulation, while response_T is the response of the NfkB system in an individual cell at time point T. The above table can be shown in R by calling

library(SLEMI)
rbind(data_nfkb[1:3,1:4],data_nfkb[10001:10003,1:4],tail(data_nfkb[,1:4],3))

2.3 Calculation of the information capacity

Calculation of the information capacity with default settings is performed by the command

capacity_logreg_main(dataRaw, signal, response, output_path)

where the required arguments are

The function returns a list with the following elements

By default, all returned elements are saved in output_path directory in a file output.rds. Along with the output data, results of the computations are visualised as the graphs listed below

2.4 Calculation of the mutual information

The function mi_logreg_main() takes a similar list of arguments and generates analogous plots to the function capacity_logreg_main(). The differences are listed below.

Firstly, user must specify the distribution of input that should be used for calculation of the mutual information. It is done by passing a numeric vector via the argument pinput of mi_logreg_main() function. Secondly, the returned list stores the value of the computed mutual information (in bits) under the element mi.

2.5 Calculation of the probabilities of correct discrimination

Calculation of the probabilities of correct discrimination between pairs of input values is performed by running the following command

prob_discr_pairwise(dataRaw, signal, response, output_path)

where the required arguments are analogous to the arguments of the functions capacity_logreg_main() and mi_logreg_main(). The probabilities of correct discrimination are computed for each pair of unique input values and returned as a list with the following elements

In addition, a plot of corresponding pie charts is created in output_path in the pdf format.

3 Diagnostic procedures

In addition to the sole calculation of the information capacity, the function capacity_logreg_main() can also be used to asses accuracy of the channel capacity estimates resulting from potentially insufficient sample size and potential over-fitting of the regression model. Two test are implemented. Precisely, the function can perform

  1. Bootstrap test - capacity is re-calculated using \(\alpha\)% of data, sampled from the original dataset without replacement. After repeating the procedure \(n\) times, standard deviation of the obtained sample can serve as an error of the capacity estimate.
  2. Over-fitting test - the original data is divided into Training and Testing datasets. Then, logistic regression is estimated using \(\alpha\)% of data (training dataset), and integrals of channel capacity are calculated via Monte Carlo using remaining \((1-\alpha)\)% of data (testing dataset). It is repeated \(n\) times.

In order to perform diagnostic tests, that by default are turned off, user must set the value of the input argument

In addition, settings of the diagnostic test can be altered by changing the following parameters

4 Additional functionalities of the function capacity_logreg_main()

In addition, to the basic functionalities described above, the function capacity_logreg_main() allows to control several other parameters of the algorithm that computes the information capacity. These parameters and their effects are listed below.

The latter two parameters, i.e lr_maxit and MaxNWts, allow to change the parameters of the logistic regression model fitting within the dependent nnet package.

5 Examples

5.1 Minimal example

Below, we present a minimal model that may serve as a quick introduction to computations within the package. Precisely, we consider a system

  1. with four different input values \(X\): 0, 0.1, 1 and 10
  2. with the conditional output, \(Y|X=x\), give by a one-dimensional log-normal distribution \(\exp\{\mathcal{N}(10\cdot\frac{x}{1+x},1)\}\)
  3. and the sample consisting of 1000 observations for each input value.

The example is analogous to the Test scenario 2 of the Supplementary Information of (Jetka et al. 2019) (Section 3.2).

Input data

Firstly, we generate a a synthetic dataset. The data corresponding to the model can be generated, and represented as the data frame tempdata with columns input and output, by running

xs=c(0,0.1,1,10) # concentration of input.
tempdata = data.frame(input = factor(c(t(replicate(1000,xs))),
                      levels=xs),
                      output =  c(matrix(rnorm(4000, mean=10*(xs/(1+xs)),sd=c(1,1,1,1)),
                                            ncol=4,byrow=TRUE) ))

The generated data.frame has the following structure

input output
1 0 1.353216
2 0 -1.052140
2001 1 3.372557
2002 1 4.941955
3999 10 9.802129
4000 10 8.544487

Calculation of the information capacity

The Information capacity can be calculated using the capacity_logreg_main() function that takes the data frame “tempdata” as dataRaw argument. Column names “input” and “output” are used as arguments signal and response, respectively. The output_path is set as “minimal_example/.” Therefore, the function is run as follows

tempoutput  <- capacity_logreg_main(dataRaw=tempdata, 
                                    signal="input", response="output", 
                                    output_path="minimal_example/")

Results of the computations are returned as a data structure described before. In addition, results are presented in the form of the following graph (by default saved as MainPlot.pdf in minimal_example/ directory). It represents the input-output data and gives the corresponding channel capacity.

Standard output graph of the minimal working example

Calculation of the mutual information

To compare mutual information of experimental data with its channel capacity, we can run (uniform distribution of input values is assumed, as default)

tempoutput_mi  <- mi_logreg_main(dataRaw=tempdata, 
                                    signal="input", response="output", 
                                    output_path="minimal_exampleMI/",
                                    pinput=rep(1/4,4)) 

and display results

print(paste("Mutual Information:", tempoutput_mi$mi,"; ",
            "Channel Capacity:", tempoutput$cc, sep=" "))
## [1] "Mutual Information: 1.48638778375912 ;  Channel Capacity: 1.54582713305531"

Alternatively, the distribution of the input can be defined with probabilities \((0.4,0.1,0.4,0.1)\)

tempoutput_mi  <- mi_logreg_main(dataRaw=tempdata, 
                                    signal="input", response="output", 
                                    output_path="minimal_exampleMI/",
                                    pinput=rc(0.4,0.1,0.4,0.1)) 

and display results

print(paste("Mutual Information:", tempoutput_mi$mi,"; ",
            "Channel Capacity:", tempoutput$cc, sep=" "))
## [1] "Mutual Information: 1.33526299432175 ;  Channel Capacity: 1.54582713305531"

Calculation of the probabilities of correct discrimination

Probabilities of correct discrimination between input values are calculated as follows

tempoutput_probs  <- prob_discr_pairwise(dataRaw=tempdata, 
                                    signal="input", response="output", 
                                    output_path="minimal_exampleProbs/") 

The above command generates graph shown in Figure 4 in the output directory

Standard output graph presenting probabilities of correct discrimination between each pair of input values.

Standard output graph presenting probabilities of correct discrimination between each pair of input values.

Diagnostics

The diagnostic test can be performed as follows

dir.create("example1_testing/")
outputCLR=capacity_logreg_main(dataRaw=data_example1,
                            signal="signal",response="response",
                            output_path="example1_testing/",
                            testing=TRUE, TestingSeed = 1234, testing_cores = 4, 
                            boot_num = 40, boot_prob = 0.8, 
                            traintest_num = 40,partition_trainfrac = 0.6)

It will run diagnostics with 40 re-sampling of the data, where bootstrap is calculated using 80% of the data, while the over-fitting test uses 60% of the original dataset.

#```{r diag1,include=FALSE}
#library(SLEMI)
#dir.create(“example1_testing/”)
#outputCLR_test=capacity_logreg_main(dataRaw=data_example1,
# signal=“signal,”response=“response,” #side_variables = NULL, testing=TRUE, TestingSeed = 1234, testing_cores = 6,
# boot_num = 30, boot_prob = 0.8,
# sidevar_num = 6,
# traintest_num = 30,partition_trainfrac = 0.6)
#```

Its results are provided in graph presented in Figure 5.

Standard output graph of the diagnostic procedures. P-values (PV) are based on empirical test either left- or right- sided. In the top axis, black dot represents the estimate of the channel capacity that involves the compete dataset, red dot is the mean of bootstrap procedures, while the bars are mean +/- sd. The remaining panels are histograms of all repetitions of a specific diagnostic procedure.

The top diagram shows the value of the capacity estimate (in black) obtained from the complete dataset and the mean value of bootstrap repetitions with indicated +/- standard deviation (in red). Plots that follow show histograms of calculated capacities for different diagnostic regimes. The black dot represents the estimate of the channel capacity based on the complete dataset. In addition, corresponding empirical p-values of both tests (left- and right-sided) are calculated to assess the randomness of obtained results (PV in the plots).

A reliable estimation of the information capacity should yield the following results of the bootstrap and overfitting tests.

  1. The bootstrap test should yield distribution of the capacity estimates with small variance. In addition, the capacity estimated based on the complete dataset should not be an outlier (p-value>0.05). Otherwise, it would indicate that the sample size is too low for an accurate estimation of the channel capacity.

  2. The over-fitting test should provide similar results. The capacity estimate obtained based on the complete dataset should lie within the distribution of capacities generated in the test. In the opposite case, it could mean that the logistic regression model does not fully grasp the essential aspects of input-output dependencies in the data.

5.2 Further step-by-step introductory examples

Two step-by-step examples that further illustrate the applicability of the SLEMI package are provided in the Section 6 of the ‘Testing procedures’ pdf file that is added to the publication (Jetka et al. 2019) and can be found here.

5.3 Examples in paper

To reproduce results of the NFkB analysis presented in the publication, see Section 7 of the ‘Testing procedures’ pdf file that is added to the publication (Jetka et al. 2019) and can be found here.

6 Session Info

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19043)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=C                           
## [2] LC_CTYPE=English_United Kingdom.1252   
## [3] LC_MONETARY=English_United Kingdom.1252
## [4] LC_NUMERIC=C                           
## [5] LC_TIME=English_United Kingdom.1252    
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
## [1] SLEMI_1.0.1    reshape2_1.4.4 stringr_1.4.0  gridExtra_2.3  ggplot2_3.3.3 
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0     xfun_0.20            purrr_0.3.4         
##  [4] splines_4.0.3        lattice_0.20-41      colorspace_2.0-0    
##  [7] vctrs_0.3.5          generics_0.1.0       stats4_4.0.3        
## [10] htmltools_0.5.1.1    yaml_2.2.1           survival_3.2-7      
## [13] prodlim_2019.11.13   rlang_0.4.9          e1071_1.7-4         
## [16] ModelMetrics_1.2.2.2 pillar_1.4.7         glue_1.4.2          
## [19] withr_2.4.0          foreach_1.5.1        lifecycle_0.2.0     
## [22] plyr_1.8.6           lava_1.6.8.1         timeDate_3043.102   
## [25] munsell_0.5.0        gtable_0.3.0         recipes_0.1.15      
## [28] codetools_0.2-18     evaluate_0.14        knitr_1.30          
## [31] caret_6.0-86         class_7.3-17         highr_0.8           
## [34] Rcpp_1.0.5           scales_1.1.1         ipred_0.9-9         
## [37] digest_0.6.27        stringi_1.5.3        dplyr_1.0.2         
## [40] tools_4.0.3          magrittr_2.0.1       tibble_3.0.4        
## [43] crayon_1.3.4         pkgconfig_2.0.3      ellipsis_0.3.1      
## [46] MASS_7.3-53          Matrix_1.2-18        data.table_1.13.4   
## [49] pROC_1.17.0.1        lubridate_1.7.9.2    gower_0.2.2         
## [52] rmarkdown_2.6        iterators_1.0.13     R6_2.5.0            
## [55] rpart_4.1-15         nnet_7.3-15          nlme_3.1-151        
## [58] compiler_4.0.3

References

Jetka, Tomasz, Karol Nienałtowski, Tomasz Winarski, Sławomir Błoński, and Michał Komorowski. 2019. “Information-Theoretic Analysis of Multivariate Single-Cell Signaling Responses.” PLOS Computational Biology 15 (7): e1007132. https://doi.org/10.1371/journal.pcbi.1007132.

  1. author and maintainer, please contact via t.jetka a t gmail.com↩︎