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
- with four different input values \(X\): 0, 0.1, 1 and 10
- 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)\}\)
- 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.
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.
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.
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.