Title: | Statistical Learning Based Estimation of Mutual Information |
---|---|
Description: | The implementation of the algorithm for estimation of mutual information and channel capacity from experimental data by classification procedures (logistic regression). Technically, it allows to estimate information-theoretic measures between finite-state input and multivariate, continuous output. Method described in Jetka et al. (2019) <doi:10.1371/journal.pcbi.1007132>. |
Authors: | Tomasz Jetka [aut, cre], Karol Nienaltowski [ctb], Michal Komorowski [ctb] |
Maintainer: | Tomasz Jetka <[email protected]> |
License: | LGPL (>=2) |
Version: | 1.0.2 |
Built: | 2024-11-13 03:28:36 UTC |
Source: | https://github.com/tjetka/slemi |
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). It is recommended to conduct estimation by calling capacity_logreg_main.R.
capacity_logreg_algorithm( data, signal = "signal", response = "response", side_variables = NULL, formula_string = NULL, model_out = TRUE, cc_maxit = 100, lr_maxit = 1000, MaxNWts = 5000 )
capacity_logreg_algorithm( data, signal = "signal", response = "response", side_variables = NULL, formula_string = NULL, model_out = TRUE, cc_maxit = 100, lr_maxit = 1000, MaxNWts = 5000 )
data |
must be a data.frame object. Cannot contain NA values. |
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 |
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. |
model_out |
is the logical indicating if the calculated logistic regression model should be included in output list |
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. |
a list with three elements:
output$cc - channel capacity in bits
output$p_opt - optimal probability distribution
output$regression - confusion matrix of logistic regression predictions
output$model - nnet object describing logistic regression model (if model_out=TRUE)
[1] 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.
tempdata=data_example1 outputCLR1=capacity_logreg_algorithm(data=tempdata, signal="signal", response="response",cc_maxit=3,model_out=FALSE, formula_string = "signal~response")
tempdata=data_example1 outputCLR1=capacity_logreg_algorithm(data=tempdata, signal="signal", response="response",cc_maxit=3,model_out=FALSE, formula_string = "signal~response")
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.
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 )
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 )
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 |
In a typical experiment aimed to quantify information flow a given signaling system, input values , ranging from 0 to saturation are considered.
Then, for each input level,
,
observations are collected, which are represented as vectors
Within information theory the degree of information transmission is measured as the mutual information
where is the marginal distribution of the output. MI is expressed in bits and
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, ,
defines the information capacity, C. Formally,
Information capacity is expressed in bits and 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,
, which is known to be a simpler problem. Estimation of the MI using estimates of
, denoted here as
, is possible as the MI, can be
alternatively written as
The expected value (as in above expression) with respect to distribution can be approximated by the average with respect to data
Here, we propose to use logistic regression as . Specifically,
Following this approach, channel capacity can be calculated by optimising MI with respect to the input distribution, .
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.
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
[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.
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
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
In the paper describing methodological aspects of our algorithm we present the analysis of information transmission
in NfkB pathway upon the stimulation of TNF-. Experimental data from this experiment in the form of single-cell
time series are attached to the package as a data.frame object and can be accessed using 'data_nfkb' variable.
Each row of ‘data_nfkb' represents a single observation of a cell. Column ’signal' indicates the level of TNF-
stimulation for a given cell, while columns 'response_T', gives the normalised ratio of nuclear and cytoplasmic transcription
factor as described in Supplementary Methods of the corresponding publication. In the CRAN version of the package
we included only a subset of the data (5 time measurements). For the full datasets, please access GitHub pages.
data_nfkb
data_nfkb
A data frame with 15632 rows and 6 variables:
Level of TNFa stimulation
The concentration of normalised NfkB transcription factor, measured at time 0
The concentration of normalised NfkB transcription factor, measured at time 3
The concentration of normalised NfkB transcription factor, measured at time 21
The concentration of normalised NfkB transcription factor, measured at time 90
The concentration of normalised NfkB transcription factor, measured at time 120
#'
For each concentration, there are at least 1000 single-cell observation (with the exception of 0.5ng stimulation, where the number of identified cells is almost 900)
in-house experimental data
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). It is recommended to conduct estimation by calling mi_logreg_main.R.
mi_logreg_algorithm( data, signal = "signal", response = "response", side_variables = NULL, pinput = NULL, formula_string = NULL, lr_maxit = 1000, MaxNWts = 5000, model_out = TRUE )
mi_logreg_algorithm( data, signal = "signal", response = "response", side_variables = NULL, pinput = NULL, formula_string = NULL, lr_maxit = 1000, MaxNWts = 5000, model_out = TRUE )
data |
must be a data.frame object. Cannot contain NA values. |
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 |
side_variables |
(optional) is a character vector that indicates side variables' columns of data, if NULL no side variables are included |
pinput |
is a numeric vector with prior probabilities of the input values. Uniform distribution is assumed as default (pinput=NULL). |
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. |
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. |
model_out |
is the logical indicating if the calculated logistic regression model should be included in output list |
a list with three elements:
output$mi - mutual information in bits
output$pinput - prior probabilities used in estimation
output$regression - confusion matrix of logistic regression model
output$model - nnet object describing logistic regression model (if model_out=TRUE)
[1] 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.
## Estimate mutual information directly temp_data=data_example1 output=mi_logreg_algorithm(data=data_example1, signal = "signal", response = "response")
## Estimate mutual information directly temp_data=data_example1 output=mi_logreg_algorithm(data=data_example1, signal = "signal", response = "response")
The main wrapping function for basic usage of SLEMI package for estimation of mutual information. 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.
mi_logreg_main( dataRaw, signal = "input", response = NULL, output_path = NULL, side_variables = NULL, pinput = NULL, formula_string = NULL, 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 )
mi_logreg_main( dataRaw, signal = "input", response = NULL, output_path = NULL, side_variables = NULL, pinput = NULL, formula_string = NULL, 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 )
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 |
pinput |
is a numeric vector with prior probabilities of the input values. Uniform distribution is assumed as default (pinput=NULL). |
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. |
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 |
In a typical experiment aimed to quantify information flow a given signaling system, input values , ranging from 0 to saturation are considered.
Then, for each input level,
,
observations are collected, which are represented as vectors
Within information theory the degree of information transmission is measured as the mutual information
where is the marginal distribution of the output. MI is expressed in bits and
can be interpreted as the number of
inputs that the system can resolve on average.
In contrast to existing approaches, instead of estimating, possibly highly dimensional, conditional output distributions , we propose to estimate the discrete, conditional input distribution,
, which is known to be a simpler problem. Estimation of the MI using estimates of
, denoted here as
, is possible as the MI, can be
alternatively written as
The expected value (as in above expression) with respect to distribution can be approximated by the average with respect to data
Here, we propose to use logistic regression as . Specifically,
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.
a list with several elements:
output$regression - confusion matrix of logistic regression predictions
output$mi - mutual information in bits
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
[1] 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.
tempdata=data_example1 outputCLR1=mi_logreg_main(dataRaw=tempdata, signal="signal", response="response") tempdata=data_example2 outputCLR2=mi_logreg_main(dataRaw=tempdata, signal="signal", response=c("X1","X2","X3")) #For further details see vignette
tempdata=data_example1 outputCLR1=mi_logreg_main(dataRaw=tempdata, signal="signal", response="response") tempdata=data_example2 outputCLR2=mi_logreg_main(dataRaw=tempdata, signal="signal", response=c("X1","X2","X3")) #For further details see vignette
Estimates probabilities of correct discrimination (PCDs) between each pair of input/signal values using a logistic regression model.
prob_discr_pairwise( dataRaw, signal = "input", response = NULL, side_variables = NULL, formula_string = NULL, output_path = NULL, scale = TRUE, lr_maxit = 1000, MaxNWts = 5000, diagnostics = TRUE )
prob_discr_pairwise( dataRaw, signal = "input", response = NULL, side_variables = NULL, formula_string = NULL, output_path = NULL, scale = TRUE, lr_maxit = 1000, MaxNWts = 5000, diagnostics = TRUE )
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 |
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. |
output_path |
is a directory where a pie chart with calculated probabilities will be saved. If NULL, the graph will not be created. |
scale |
is a logical indicating if the response variables should be scaled and centered before fitting logistic regression |
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. |
diagnostics |
is a logical indicating if details of logistic regression fitting should be included in output list |
In order to estimate PCDs, for a given pair of input values and
, we propose to fit a logistic regression model
using response data corresponding to the two considered inputs, i.e.
, for
and
ranging from
1 to
.
To ensure that both inputs have equal contribution to the calculated discriminability, equal probabilities should be assigned,
. Once the regression model is fitted, probability of assigning a given cellular response,
,
to the correct input value is estimated as
Note that as well as
The average of the above probabilities over all observations
yields PCDs
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.
a list with two elements:
output$prob_matr - a matrix, where
is the number of inputs, with probabilities of correct
discrimination between pairs of input values.
output$diagnostics - (if diagnostics=TRUE) a list corresponding to logistic regression models fitted for each pair of input values. Each element consists of three sub-elements: 1) nnet_model - nnet object summarising logistic regression model; 2) prob_lr - probabilities of assignment obtained from logistic regression model; 3) confusion_matrix - confusion matrix of classificator.
[1] 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.
## Calculate probabilities of discrimination for nfkb dataset it=21 # choose from 0, 3, 6, ..., 120 for measurements at other time points output=prob_discr_pairwise(dataRaw=data_nfkb[data_nfkb$signal%in%c("0ng","1ng","100ng"),], signal = "signal", response = paste0("response_",it))
## Calculate probabilities of discrimination for nfkb dataset it=21 # choose from 0, 3, 6, ..., 120 for measurements at other time points output=prob_discr_pairwise(dataRaw=data_nfkb[data_nfkb$signal%in%c("0ng","1ng","100ng"),], signal = "signal", response = paste0("response_",it))