# Algorithm-agnostic significance testing in supervised learning with
multimodal data

Lucas Kook∗, Anton Rask Lundborg

(Department of Mathematical Sciences, University of Copenhagen, Denmark
∗E-mail<EMAIL_ADDRESS>
###### Abstract

Valid statistical inference is crucial for decision-making but difficult to
obtain in supervised learning with multimodal data, e.g., combinations of
clinical features, genomic data, and medical images. Multimodal data often
warrants the use of black-box algorithms, for instance, random forests or
neural networks, which impede the use of traditional variable significance
tests. We address this problem by proposing the use of COvariance Measure
Tests (COMETs), which are calibrated and powerful tests that can be combined
with any sufficiently predictive supervised learning algorithm. We apply
COMETs to several high-dimensional, multimodal data sets to illustrate (i)
variable significance testing for finding relevant mutations modulating drug-
activity, (ii) modality selection for predicting survival in liver cancer
patients with multiomics data, and (iii) modality selection with clinical
features and medical imaging data. In all applications, COMETs yield results
consistent with domain knowledge without requiring data-driven pre-processing
which may invalidate type I error control. These novel applications with high-
dimensional multimodal data corroborate prior results on the power and
robustness of COMETs for significance testing. The comets R package and source
code for reproducing all results is available at
https://github.com/LucasKook/comets. All data sets used in this work are
openly available.

##  1 Introduction

A fundamental challenge of modern bioinformatics is dealing with the
increasingly multimodal nature of data (Cheerla and Gevaert, 2019; Ahmed et
al., 2021; Stahlschmidt et al., 2022). The task of _supervised learning_ ,
that is, the problem of predicting a response variable $Y$ from features $X$,
has received considerable attention in recent years resulting in a plethora of
algorithms for a wide range of settings that permit prediction using several
data modalities simultaneously (Hastie et al., 2009). With the advent of deep
learning, even non-tabular data modalities, such as text or image data, can be
included without requiring manual feature engineering (LeCun et al., 2015).
Methods such as these are highly regularized (if trained correctly) which
minimizes the statistical price of adding too many irrelevant variables.
However, continuing to collect features or modalities that do not contribute
to the predictiveness of a model still has an economic cost and, perhaps more
importantly, it is of scientific interest to determine whether a particular
feature or modality $X$ adds predictive power in the presence of additional
features or modalities $Z$ (Smucler and Rotnitzky, 2022).

The problem of determining which features or modalities are significantly
associated with the response is usually addressed by means of _conditional
independence testing_. The response $Y$ is independent of the modality $X$
given further modalities $Z$ if the probability that $Y$ takes any particular
value knowing both $X$ and $Z$ is the same as the probability knowing just
$Z$. In particular, $X$ does not help in predicting $Y$ if $Z$ is taken into
account already (see Section 2.1 for a more precise definition).

Traditional variable significance tests start by posing a parametric
relationship between the response $Y$ and features $X$ and $Z$, for instance,
the Wald test in a generalized linear model. When $X$ or $Z$ are complicated
data modalities, it is seldom possible to write down a realistic model for
their relationship with $Y$; thus a different approach is required.
Furthermore, even when models can be explicitly parametrized, it is not clear
that the resulting tests remain valid when the model is not specified
correctly (Shah and Bühlmann, 2023).

More recently, kernel-based conditional independence tests have been proposed
which use a characterization of conditional independence by means of kernel
embeddings to construct tests (Zhang et al., 2012; Strobl et al., 2019).
However, these tests are difficult to calibrate in practice and rely
intimately on kernel ridge regression. Several alternative algorithm-agnostic
tests have been developed under the so-called ‘Model-X’ assumption where one
supposes that a model is known (or at least estimable to high precision) for
the full distribution of $X$ given $Z$ (Candès et al., 2018; Berrett et al.,
2019). Given the difficulty of learning conditional distributions such an
assumption is rarely tenable. Algorithm-agnostic variable importance measures
have also been developed with statistically optimal estimators (Williamson et
al., 2021, 2023). However, efficient estimation of an importance measure does
not necessarily translate to an optimal test to distinguish between
conditional dependence and independence (see, e.g., the introduction of
Lundborg et al., 2022a).

Figure 1: Overview of the proposed algorithm-agnostic significance testing
framework for multimodal data using COMETs. Variable significance:
Differential gene expression can be assessed in presence of the potentially
high-dimensional/non-tabular confounder $Z$. Modality selection: Entire
modalities can be subjected to significance testing, which lends itself to
modality selection in multi-omics applications. MSPE: Mean squared prediction
error.

In this paper, we describe a family of significance tests referred to
collectively as COvariance Measure Tests (COMETs) which are algorithm-agnostic
and valid (in the sense of controlling the probability of false positives) as
long as the algorithms employed are sufficiently predictive. We will primarily
focus on the Generalised Covariance Measure (GCM) test (Shah and Peters, 2020)
which we think of as an ‘all-purpose’ test that should be well-behaved in most
scenarios and the more complicated Projected Covariance Measure (PCM) test
which is more flexible but may require a more careful choice of algorithms.
Figure 1 gives an overview of the proposed algorithm-agnostic significance
testing framework based on COMETs and the types of applications that are
presented in this manuscript.

##  2 Methods

In this section, we first provide some background on conditional independence.
We then move on to describe the computation of the GCM and PCM tests in
addition to the assumptions required for their validity. Finally, we describe
the datasets which we will analyze in Section 3.

###  2.1 Background on conditional independence

For a real-valued response $Y$ and features $X$ and $Z$, we say that $Y$ is
_conditionally independent of $X$ given $Z$_ and write
$Y\perp\nolinebreak\perp X\,|\,Z$ if

$\mathbb{E}[h(Y)\,|\,X,Z]=\mathbb{E}[h(Y)\,|\,Z]\quad\text{for all functions
$h$.}$ (1)

That is, for any transformation $h$ of $Y$, the best predictor (in a mean-
squared error sense) of $h(Y)$ using both $X$ and $Z$ is equal to the best
predictor using just $Z$.111An alternative characterization, when $Y$ has a
conditional density given $X$ and $Z$ denoted by $f_{Y|X,Z}(y|x,z)$, is given
by: $f_{Y|X,Z}(y|x,z)=f_{Y|Z}(y|z)$ if and only if $Y$ is independent of $X$
given $Z$.

A helpful starting point for the construction of a conditional independence
test is to consider the product of a population residual from a $Y$ on $Z$
regression $\varepsilon:=Y-\mathbb{E}[Y\,|\,Z]$ and, for now considering a
one-dimensional $X$, from an $X$ on $Z$ regression
$\xi:=X-\mathbb{E}[X\,|\,Z]$. As these are population residuals, $Z$ is no
longer helpful in predicting their values, so
$\mathbb{E}[\varepsilon\,|\,Z]=\mathbb{E}[\varepsilon]=0$ and similarly
$\mathbb{E}[\xi\,|\,Z]=0$. When $Y\perp\nolinebreak\perp X\,|\,Z$, we can say
more: the product of the residuals is also mean zero since

$\displaystyle\mathbb{E}[\varepsilon\xi]=\mathbb{E}[\mathbb{E}[\varepsilon\xi\,|\,X,Z]]=\mathbb{E}[\mathbb{E}[\varepsilon\,|\,X,Z]\xi]=\mathbb{E}[\mathbb{E}[\varepsilon\,|\,Z]\xi]=0,$
(2)

where the second equality uses that $\xi$ is perfectly predicted using $X$ and
$Z$, the third equality uses (1) with $h(y)=y$ and the final equality uses
$\mathbb{E}[\varepsilon\,|\,Z]=0$. The GCM test is based on testing whether
$\mathbb{E}[\varepsilon\xi]=0$ and we will describe the details of how to
compute it in Section 2.2. For the GCM test to perform well, it is important
to determine when we can expect $\mathbb{E}[\varepsilon\xi]$ to be non-zero
under conditional dependence. When $Y$ follows a partially linear model given
$X$ and $Z$, that is, $\mathbb{E}[Y\,|\,X,Z]=\theta X+h(Z)$ for some function
$h$, then $\mathbb{E}[\varepsilon\xi]\neq 0$ exactly when $\theta\neq 0$ and
the magnitude of $\mathbb{E}[\varepsilon\xi]$ is proportional to $\theta$.
This includes as a special case the linear model for $Y$ given $X$ and $Z$.
There is a natural generalization of (2) to the case where $X$ is a
$d$-dimensional vector, where the equation is interpreted component-wise in
$X$. Although the GCM is also defined in these settings, computing the test
involves many regressions when $X$ is high-dimensional which can be
impractical (see Section 2.2.3).

Unfortunately, it is not difficult to come up with examples where
$Y\not\perp\nolinebreak\perp X\,|\,Z$ but $\mathbb{E}[\varepsilon\xi]=0$. For
instance, if $X$ and $Z$ are independent and standard normally distributed and
$\mathbb{E}[Y\,|\,X,Z]=X^{2}$, then $\mathbb{E}[Y\,|\,Z]=1$ (since $Z$ carries
no information about $X$ so the best predictor is just the mean of $Y$) hence

$\mathbb{E}[\varepsilon\xi]=\mathbb{E}[(X^{2}-1)X]=0,$

using that $\mathbb{E}[X]=\mathbb{E}[X^{3}]=0$ for a standard normal variable.
A more elaborate example is given in Figure 2 (left and middle panel) and even
more examples exist when $X$ and $Z$ are dependent (see Lundborg et al. 2022a,
Section 6 and Scheidegger et al. 2022, Section 3.1.2). We now describe a test
which can detect such dependencies.

Figure 2: Illustration of the GCM and PCM test under the alternative that
$Y\in\mathbb{R}$ is not conditionally independent of $X\in\mathbb{R}$ given
$Z\in\mathbb{R}$, where $Y\coloneqq f(X)g(Z)+\varepsilon$ and
$f(x):=1+\sin(3x^{2})$, $g(z):=1+z^{3}$. The GCM test first computes the
residual for the regression of $Y$ on $Z$, which shows no correlation with
$X$. Thus, the GCM fails to reject (correlation coefficient $R=0.01$). The
PCM, in addition, learns the optimal transformation of $X$ (depending on $Z$)
to test conditional mean independence of $Y$ and $X$ given $Z$. Thus, in this
example, the PCM test correctly rejects ($R=0.87$). Although the residuals in
the second panel are clearly not independent, it is not valid to conclude
conditional dependence from rejecting an independence test here (see Shah and
Peters, 2020, Example 1).

A more ambitious target is to detect whenever an arbitrary (e.g., non-tabular)
$X$ is helpful for the prediction of $Y$ in the presence of $Z$ measured in
terms of mean-squared error. To achieve this goal, we can use the fact,
derived in the same way as (2), that

$\mathbb{E}[\varepsilon f(X,Z)]=0\quad\text{for all functions $f$},$

whenever $Y\perp\nolinebreak\perp X\,|\,Z$. The GCM targets the quantity
involving the function $f(X,Z)=X$. However, by instead using
$f(X,Z)=\mathbb{E}[Y\,|\,X,Z]-\mathbb{E}[Y\,|\,Z]$ (which depends on the joint
distribution of $X$ and $Z$), we obtain that

$\mathbb{E}[\varepsilon
f(X,Z)]=\mathbb{E}[(\mathbb{E}[Y\,|\,Z]-\mathbb{E}[Y\,|\,X,Z])^{2}]=:\tau.$
(3)

This quantity is strictly greater than $0$ if and only if $X$ is helpful for
the prediction of $Y$ in the presence of $Z$. The PCM test is based on testing
whether $\tau=0$ and we will describe the details of how to compute it in
Section 2.2. In fact, the PCM is based on an alternative $f$ given by
$f(X,Z)=(\mathbb{E}[Y\,|\,X,Z]-\mathbb{E}[Y\,|\,Z])/\operatorname{Var}(Y\,|\,X,Z)$
that turns out to result in a more powerful test (see Figure 2 and Lundborg et
al., 2022a, Section 1.1). An added benefit of tests targeting $\tau$ is that
no regressions are needed with $X$ as the response, which can vastly reduce
the computational burden when compared to tests that target
$\mathbb{E}[\varepsilon\xi]$.

The targets mentioned above rely intimately on population quantities that are
unknown and hence need to be estimated when computing tests in practice. To
ensure that the estimation errors do not interfere with the performance of the
tests, we need to be able to learn the functions to a sufficient degree of
accuracy. These requirements put restrictions on when the GCM and PCM are
valid tests but such restrictions are not unique to these tests. In fact,
unless $Z$ is discrete, it is impossible to construct an assumption-free
conditional independence test that simultaneously controls the probability of
false rejections and is able to detect dependence (Shah and Peters, 2020; Kim
et al., 2022). This result implies that additional assumptions need to be
imposed to ensure the feasibility of testing for conditional independence.

###  2.2 Covariance measure tests

We now describe the specifics of computing the GCM and the PCM. For the
remainder of this section, we assume that we have a dataset consisting of $n$
independent observations of a real-valued response $Y$ and some additional
features or modalities $X$ and $Z$.

####  2.2.1 Generalised Covariance Measure (GCM) test

The GCM test is based on (2) but to compute the test in practice, we need to
form an empirical version of the equation. For simplicity, we consider, for
now, $X\in\mathbb{R}$. Let $\widehat{\xi}_{i}$ denote the residual for the
$i$th observation from regressing $X$ on $Z$ and similarly
$\widehat{\varepsilon}_{i}$ from regressing $Y$ on $Z$. We now test
$H_{0}:\mathbb{E}[\varepsilon\xi]=0$ by comparing

$T:=\frac{\left(\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}\widehat{\xi}_{i}\right)^{2}}{\frac{1}{n}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}^{2}\widehat{\xi}_{i}^{2}-\left(\frac{1}{n}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}\widehat{\xi}_{i}\right)^{2}}$
(4)

to a $\chi^{2}_{1}$ distribution. The term inside the square in the numerator
is $\sqrt{n}$ times an estimate of (2) while the denominator standardizes the
variance of the test statistic. The test statistic in (4) is approximately
$\chi^{2}_{1}$ for large enough sample sizes if the regression methods
employed are sufficiently predictive and $Y\perp\nolinebreak\perp X\,|\,Z$
(Shah and Peters, 2020, Theorem 6). Note that the procedure above did not use
anything special about $Z$ other than the existence of a regression method
that can approximate the conditional expectations of $Y$ and $X$ given $Z$.
The computations above naturally generalize to settings where
$X\in\mathbb{R}^{d}$ and we summarize the general procedure in Algorithm 1.

Algorithm 1 GCM test of $H_{0}:\mathbb{E}[\varepsilon\xi]=0$

1:Data set containing $n$ realizations of $Y\in\mathbb{R}$,
$X\in\mathbb{R}^{d}$, $Z\in\mathcal{Z}$, algorithms for the regressions of $Y$
on $Z$ and $X_{j}$ on $Z$ for $j\in\\{1,\dots,d\\}$

2:Run the regression of $Y$ on $Z$ and compute residuals
$\widehat{\varepsilon}_{i}$

3:Run the regressions of $X_{j}$ on $Z$ for $j\in\\{1,\dots,d\\}$
and compute residuals
$\widehat{\xi}_{i}\leftarrow(\widehat{\xi}_{i1},\dots,\widehat{\xi}_{id})$

4:Compute the mean residual product:
$L\leftarrow\frac{1}{n}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}\widehat{\xi}_{i}$

5:Compute the residual variance:
$\Sigma\leftarrow\frac{1}{n}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}^{2}\widehat{\xi}_{i}\widehat{\xi}_{i}^{\top}-LL^{\top}$

6:Compute the test statistic: $T\leftarrow
n\lVert\Sigma^{-1/2}L\rVert_{2}^{2}$

7:Compute the $p$-value: $p\leftarrow 1-F_{\chi^{2}_{d}}(T)$

8:return $p$

####  2.2.2 Projected Covariance Measure (PCM) test

The computation of the PCM test is more challenging than the computation of
the GCM test since the PCM requires learning
$f(X,Z)=(\mathbb{E}[Y\,|\,X,Z]-\mathbb{E}[Y\,|\,Z])/\operatorname{Var}(Y\,|\,X,Z)$
to be able to estimate $\mathbb{E}[\varepsilon f(X,Z)]$. Furthermore, $f$
cannot be learned on the same observations that are used to compute the test
statistic as this would potentially result in dependence between the residuals
constituting the test statistic and thus in many false rejections when
$Y\perp\nolinebreak\perp X\,|\,Z$.

The first step when computing the test statistic of the PCM test is therefore
to split the dataset in two halves $D_{1}$ and $D_{2}$ of equal size (for
simplicity, we assume that we have $2n$ observations, so both $D_{1}$ and
$D_{2}$ are of size $n$). On $D_{2}$, we compute an estimate $\widehat{f}$ of
$f$ by first regressing $Y$ on $X$ and $Z$ yielding an estimate $\widehat{g}$
and regressing $Y$ on $Z$ yielding an estimate $\widehat{m}$. We then regress
$(Y-\widehat{g}(X,Z))^{2}$ on $(X,Z)$ on $D_{2}$ yielding an estimate of
$\operatorname{Var}(Y\,|\,X,Z)$ which we denote $\widehat{v}$. We now set
$\widehat{f}(x,z):=(\widehat{g}(x,z)-\widehat{m}(z))/\widehat{v}(x,z)$ and,
working on $D_{1}$, we regress $Y$ on $Z$ yielding a residual for the $i$th
observation $\widehat{\varepsilon}_{i}$ and we regress $\widehat{f}(X,Z)$ on
$Z$ yielding a residual $\widehat{\zeta}_{i}$. Finally, we compute

$T:=\frac{\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}\widehat{\zeta}_{i}}{\left(\frac{1}{n}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}^{2}\widehat{\zeta}_{i}^{2}-\left(\frac{1}{n}\sum_{i=1}^{n}\widehat{\varepsilon}_{i}\widehat{\zeta}_{i}\right)^{2}\right)^{1/2}}$
(5)

and reject the null by comparing to a standard normal distribution. In fact,
as the target of $T$ in (3) is positive under conditional dependence, we
perform a one-sided test which rejects when $T$ is large. The test statistic
in (5) is approximately standard Gaussian if the regression methods employed
for the $\widehat{f}(X,Z)$ on $Z$ and $Y$ on $Z$ are sufficiently predictive,
the estimates $\widehat{f}$ are not too complicated and
$Y\perp\nolinebreak\perp X\,|\,Z$ (Lundborg et al., 2022a, Theorem 4). The
test is powerful against alternatives where $\widehat{f}$ is correlated with
the true $f$ and the aforementioned regression methods remain powerful
(Lundborg et al., 2022a, Theorem 5). We summarize the procedure in Algorithm 2
below.222In this description and in Algorithm 2 we have omitted a few minor
corrections to the estimation of $\widehat{f}$ that are done for numerical
stability or as finite sample corrections. The full version of the algorithm
with these additions is given in Lundborg et al. (2022a, Algorithm 1).

Algorithm 2 PCM test of $H_{0}:\mathbb{E}[Y\,|\,X,Z]=\mathbb{E}[Y\,|\,Z]$

1:Data containing $2n$ realizations of $Y\in\mathbb{R}$, $X\in\mathcal{X}$,
$Z\in\mathcal{Z}$, regression algorithms as outlined in Section 2.2.2

2:Split the data randomly into equal parts $D_{1}$ and $D_{2}$ of size $n$

3:Regress $Y$ on $X,Z$ using $D_{2}$ to obtain $\widehat{g}$

4:Regress of $Y$ on $Z$ using $D_{2}$ to obtain $\widehat{m}$

5:Regress $(Y-\widehat{g}(X,Z))^{2}$ on $X,Z$ using $D_{2}$ to obtain
$\widehat{v}$

6:Set
$\widehat{f}(x,z)\leftarrow(\widehat{g}(x,z)-\widehat{m}(z))/\widehat{v}(x,z)$

7:Regress $Y$ on $Z$ using $D_{1}$ to compute residuals
$\widehat{\varepsilon}_{i}$

8:Regress $\widehat{f}(X,Z)$ on $Z$ using $D_{1}$ to compute residuals
$\widehat{\zeta}_{i}$

9:Compute $T$ as in (5) using $\widehat{\varepsilon}_{i}$ and
$\widehat{\zeta}_{i}$, $i\in\\{1,\dots,n\\}$

10:Compute $p$-value: $p\leftarrow 1-\Phi(T)$

11:return $p$

Due to the sample splitting, the $p$-value of the PCM is a random quantity. We
can compute the PCM on several different splits to produce multiple $p$-values
that can be dealt with using standard corrections for multiple testing. In
practice, we instead compute the $p$-value as in step 9 of Algorithm 2 but
instead using the average of the test statistics from the different splits. We
denote the number of different splits by $K$ and use 5–10 in the applications.
The resulting test should be conservative which results in a power loss,
however, the test averaged from different splits should still be more powerful
than a single application of the PCM due to more efficient use of the data. If
one desires a perfectly calibrated $p$-value from multiple splits, it is
possible to use the method in Guo and Shah (2023) but we do not pursue this
further here.

####  2.2.3 Comparison of the GCM and PCM tests

The GCM and PCM tests not only differ in terms of their target quantities, but
also regarding computational aspects. The GCM test requires the regression of
$Y$ on $Z$ and $X$ on $Z$. This prohibits the use of the GCM in settings where
$X$ is a high-dimensional or non-tabular data modality and can not be
represented as or reduced to a low-dimensional tabular modality. The PCM test,
on the other hand, does not require regressing $X$ on $Z$. Thus, the PCM test
allows the end-to-end use of non-tabular data modalities, such as images or
text, for instance, via the use of deep neural networks. In contrast to the
GCM, the PCM relies on sample splitting and requires more regressions and may
thus be less data-efficient. This is addressed, in parts, by repeating the PCM
test with multiple random splits, as described above.

###  2.3 Data sets

####  2.3.1 Variable significance testing: CCLE data

We consider a subset of the anti-cancer drug dataset from the Cancer Cell Line
Encyclopedia (CCLE, Barretina et al., 2012) which contains the response to the
PLX4720 drug as a summary measure obtained from a dose-response curve and a
set of $1638$ mutations (absence/presence) in $n=472$ cancer cell lines.
Following the pre-processing steps in Bellot and van der Schaar (2019) and Shi
et al. (2021), we screen for mutations which are marginally correlated with
drug response
$S\coloneqq\\{j\in[1638]:\lvert\operatorname{Cor}(Y,X_{j})\rvert>0.05\\}$,
which leaves $\lvert S\rvert=466$ mutations.

####  2.3.2 Modality selection: TCGA data

We consider the openly available TCGA HCC multiomics data set used in
Chaudhary et al. (2018); Poirion et al. (2021). The preprocessed data consist
of survival times for $n=360$ patients with liver cancer together with RNA-seq
($\operatorname{RNAseq}\in\mathbb{R}^{15629}$), miRNA
($\operatorname{miRNA}\in\mathbb{R}^{365}$), and DNA methylation
($\operatorname{DNAm}\in\mathbb{R}^{19883}$) modalities. Details on the pre-
processing steps can be found in Chaudhary et al. (2018).

####  2.3.3 Modality selection with imaging: MIMIC data

We consider the MIMIC Chest X-Ray data set (Johnson et al., 2019; Sellergren
et al., 2022), which contains the race ($R$; with levels “white”, “black”,
“asian”), sex ($X_{1}$; with levels “male”, “female”), age ($X_{2}$, in
years), pre-trained embeddings of chest x-rays ($Z$) and (among other response
variables) whether a pleural effusion ($Y$) was visible on the x-ray for
$n=181342$ patients. The dimension of the image embedding was reduced by using
the first 111 components of a singular value decomposition, which explain 98%
of the variance.

###  2.4 Computational details

All analyses were carried out using the R language for statistical computing
(R Core Team, 2021). The covariance measure tests are implemented in comets
which relies on ranger (Wright and Ziegler, 2017) and glmnet (Tay et al.,
2023) for the random forest and LASSO regressions, respectively. Code for
reproducing all results is available at https://github.com/LucasKook/comets.
In the following, unless specified otherwise, GCM and PCM tests are run with
random forests for all regressions.

Table 1: Results for the CCLE data in Section 3.1. The table shows variable
importance ranks and $p$-values for the relation of mutations of ten genes
with the response to PLX4720 conditional on the 465 other mutations in the
data. The PCM test was run with $K=10$ random splits. The variable importance
ranks and the CRT, GCIT and DGCIT results were obtained from Bellot and van
der Schaar (2019) and Shi et al. (2021).

Method | Gene mutations
---|---
| BRAF_V600E | BRAF_MC | HIP1 | FLT3 | CDC42BPA | THBS3 | DNMT1 | PRKD1 | PIP5K1A | MAP3K5
EN | 1 | 3 | 4 | 5 | 7 | 8 | 9 | 10 | 19 | 78
RF | 1 | 2 | 3 | 14 | 8 | 34 | 28 | 18 | 7 | 9
CRT | $<0.001$ | $<0.001$ | $0.008$ | 0.017 | 0.009 | 0.017 | 0.022 | 0.002 | 0.024 | 0.012
GCIT | $<0.001$ | $<0.001$ | $0.008$ | 0.521 | 0.050 | 0.013 | 0.020 | 0.002 | 0.001 | $<0.001$
DGCIT | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | $<0.001$ | 0.794
GCM | 0.030 | 0.033 | 0.010 | 0.005 | 0.004 | 0.042 | 0.010 | 0.165 | 0.464 | 0.504
PCM | 0.001 | 0.012 | 0.008 | 0.009 | 0.014 | 0.027 | 0.014 | 0.011 | 0.022 | 0.022
GCM (no screening) | 0.003 | 0.021 | 0.006 | 0.002 | 0.002 | 0.068 | 0.007 | 0.007 | 0.223 | 0.216
PCM (no screening) | 0.002 | 0.007 | 0.082 | 0.151 | 0.186 | 0.134 | 0.138 | 0.108 | 0.198 | 0.122

##  3 Results

With our analyses, we aim to show how testing with covariance measures can be
used to tackle two of the most common supervised learning problems in
biomedical applications with multimodal data: Variable significance testing
and modality selection (see Figure 1). Throughout, we compare COMETs with
existing methods (if applicable) on openly available real data sets (see
Section 2.3 for an overview of the data sets).

###  3.1 Variable significance testing

We apply COMETs to the anti-cancer drug dataset from the Cancer Cell Line
Encyclopedia (Barretina et al., 2012) and compare with the results obtained
using the CRT (Candès et al., 2018) GCIT (Bellot and van der Schaar, 2019) and
DGCIT (Shi et al., 2021). The null hypotheses
$H_{0}(j):Y\perp\nolinebreak\perp X_{j}\,|\,X_{-j}$ are tested for $j\in S$ to
detect mutations that are significantly associated with PLX4720 drug response.

####  3.1.1 COMETs identify mutations associated with PLX4720 drug activity

Table 1 summarizes the results for the GCIT, DGCIT, GCM, and PCM test and the
ten selected mutations in Bellot and van der Schaar (2019, Figure 4). Overall,
there is large agreement between all tests which all reject the null
hypothesis for the BRAF_V600E, BRAF_MC, HIP1, FLT3, THBS3 and DNMT1 mutations,
corroborating previously reported results. For the PRKD1, PIP5K1A and MAP3K5
mutations, the PCM test rejects while the GCM test does not, which is
consistent with the PCM test having power against a larger class of
alternatives (Figure 2).

####  3.1.2 COMETs detect relevant mutations without pre-screening

Prior results rely on pre-screening genes based on their marginal correlation
with the drug response. However, marginal correlation cannot inform subsequent
conditional independence tests in general and the data-driven pre-screening
may have lead to inflated false positive rates (Berk et al., 2013). However,
the GCM and PCM test can be applied without pre-screening and still
consistently reject the null hypothesis of conditional independence for the
BRAF_V600E and BRAF_MC mutations (see rows in Table 1 with “no screening”).
When correcting (Holm) the $p$-values to attain a family-wise error rate of
$5\%$ for the ten mutations of interest, the GCM and PCM still reject the null
hypothesis for BRAF_V600E ($p=0.024$ for the GCM test and $p=0.020$ for the
PCM test). This rejection is expected because PLX4720 was designed as a BRAF
inhibitor (Barretina et al., 2012).

###  3.2 Modality selection

The goal of our analysis is to identify modalities among RNA-seq, miRNA and
DNA methylation that are important for predicting survival of liver cancer
patients by testing if the event is independent of the modality $M_{j}$ given
the other modalities $M_{-j}$, $j\in\\{1,2,3\\}$. This is a challenging
problem due to the high dimensionality of both the candidate modality $M_{j}$
and the conditioning variables in $M_{-j}$.

####  3.2.1 Evidence of DNA methylation being important for predicting
survival in liver cancer patients

Table 2 (PCM-RF) shows $p$-values for the PCM test ($K=10$ different splits)
testing for significance of the RNA-seq, miRNA and DNA methylation modalities
conditional on the remaining two without pre-screening features in any of the
modalities using a random forest regression. There is some evidence that the
DNA methylation modality is important for predicting death in liver cancer
patients. Conversely, the PCM test does not provide evidence that survival
depends on the RNA-seq or miRNA modalities, when already conditioning on the
DNA methylation data. Comparable results are obtained when substituting the
random forest regression for $Y$ on $\operatorname{RNAseq}$,
$\operatorname{miRNA}$, and $\operatorname{DNAm}$ with a cross-validated LASSO
regression using the optimal tuning parameter: After a multiple testing
correction (Holm), both PCM tests reject the null hypothesis only for the DNA
methylation modality.

Table 2: Results ($p$-values) for the multiomics application in Section 3.2
using the PCM with $K=10$ random splits once using a random forest (RF) for
the regression of $Y$ on $\operatorname{RNAseq}$, $\operatorname{miRNA}$ and
$\operatorname{DNAm}$, and once a cross-validated high-dimensional linear
regression (LASSO).

Null hypothesis | PCM-RF | PCM-LASSO
---|---|---
$H_{0}:Y\perp\nolinebreak\perp\operatorname{RNAseq}\,|\,\operatorname{miRNA},\operatorname{DNAm}$ | 0.178 | 0.066
$H_{0}:Y\perp\nolinebreak\perp\operatorname{miRNA}\,|\,\operatorname{RNAseq},\operatorname{DNAm}$ | 0.165 | 0.044
$H_{0}:Y\perp\nolinebreak\perp\operatorname{DNAm}\,|\,\operatorname{RNAseq},\operatorname{miRNA}$ | 0.014 | 0.002

###  3.3 Modality selection with imaging data

Glocker et al. (2023) provide evidence that both the race and the response
(pleural effusion) can be predicted from the x-ray embedding with high
accuracy. The goal of our analysis is to test, whether race helps predict the
response when already conditioning on age, sex, and the x-rays and, vice
versa, whether the x-rays contain information for predicting pleural effusion
given sex, age, and race.

####  3.3.1 Strong evidence for x-ray imaging and race being important for
predicting pleural effusion

There is strong evidence against the null hypotheses of pleural effusion being
independent of either x-ray imaging or race given the other and, additionally,
sex and age of a patient (Table 3).

Table 3: Results ($-\log_{10}$-transformed $p$-values) for the GCM and PCM applied to the full MIMIC data set in Section 3.3. Both tests reject both hypotheses. See Figure 3 for an uncertainty assessment.  Null hypothesis | GCM | PCM
---|---|---
$H_{0}:Y\perp\nolinebreak\perp R\,|\,X,Z$ | 6.158 | 77.762
$H_{0}:Y\perp\nolinebreak\perp Z\,|\,X,R$ | 13805.802 | 1270.361

To gauge the uncertainty in the results of the COMETs, we repeat the tests on
75 random (non-overlapping) subsamples of different sample sizes (150, 600,
2400) of the data. Only the PCM rejects the null hypothesis of pleural
effusion (PE) being independent of race given the x-ray, sex, and age of a
patient at any of the considered sample sizes which provides evidence that
$\mathbb{E}[\varepsilon\xi]$ is close to zero yet $\mathbb{E}[Y\,|\,X,Z]$
still varies non-linearly with $X$. At full sample size, the GCM does reject,
indicating the presence of a weak linear signal (estimated correlations
between pleural effusion and race residuals are smaller than $0.015$). It is
somewhat unsurprising to see both COMETs reject the null hypothesis at such
large sample sizes (Greenland, 2019, Section 4.3).

Both tests reject the null hypothesis of pleural effusion (PE) being
independent of x-ray given race, sex and age of a patient at any sample size
but in fact the GCM produces smaller $p$-values. This indicates that there is
a significant component in $\mathbb{E}[Y\,|\,X,Z]$ varying linearly with $X$;
in these cases the PCM will not outperform the GCM for a fixed sample size.

Figure 3:  Results ($-\log_{10}$-transformed $p$-values) for the GCM and PCM
applied to 75 random non-overlapping splits of different sample sizes
($n\in\\{150,600,2400\\}$) of the MIMIC data set in Section 3.3. Splitting the
data enables an analysis of the uncertainty in the tests’ rejections and the
strength of evidence against the null.

##  4 Discussion

We present COMETs for algorithm-agnostic significance testing with multi-
modal, potentially non-tabular data, which relies on tests of conditional
independence based on covariance measures. The versatility of the GCM and PCM
tests is shown in several applications involving variable significance testing
and modality selection in the presence of high-dimensional conditioning
variables. In the following, we discuss the applications in more detail and
end with a discussion of computational aspects and recommendations for using
COMETs in supervised learning applications with multimodal data.

###  4.1 Variable significance testing

The GCM and PCM test show comparable results to competing methods and can be
applied without relying on data-driven pre-screening which otherwise can
invalidate $p$-values and lead to inflated type I error rates. Type I error
control additionally suffers from the performed number of tests. After
correcting for multiple testing, the COMETs provide evidence that BRAF_V600E
is associated with PLX4720 activity while controlling for all other mutations.
As highlighted before, this is expected since PLX4720 was designed as a BRAF
inhibitor.

###  4.2 Modality selection

The PCM test is applied to the TGCA data set to test which modalities (RNAseq,
miRNA, DNAm) are important (conditional on the others) for predicting survival
in liver cancer patients and rejects the null hypothesis for the DNA
methylation modality. Failure to reject the null hypothesis for the RNA-seq
and miRNA modalities may be due to the low sample size and extremely high
dimensionality of the problem and ought to be interpreted as lack of evidence
that RNA-seq and miRNA data contain information for predicting survival beyond
DNA methylation in the data at hand. Taken together, this application
demonstrates that the PCM test can be used for modality selection with high-
dimensional candidate and conditioning modalities. COMETs could, for instance,
be used to trade off the economic cost of measuring an omics (or imaging, as
in the MIMIC application) modality with the gain in predictive power at a
given significance level.

###  4.3 Modality selection with imaging data

The large and openly available MIMIC data set serves as an example application
of how image and other non-tabular modalities may enter an analysis based on
covariance measure tests. The PCM does not require pre-trained embeddings and
could, in principle, also be used in combination with deep convolutional
neural networks if the raw imaging data is available. The 111-dimensional
embedding further enables the use of the GCM test to serve as a benchmark.
However, it is important to properly choose the regressions involved in COMETs
as the tests rely on their quality and asymptotic properties (Shah and Peters,
2020; Lundborg et al., 2022a). Nevertheless, to the best our knowledge, no
other tests exist with theoretical guarantees that also permit testing
$Y\perp\nolinebreak\perp X\,|\,Z$ when $X$ is a non-tabular modality.

###  4.4 Recommendations and outlook

As outlined in Section 2.2.3, the regression of $X$ on $Z$ required by the GCM
can become computationally challenging if $X$ is high-dimensional (which is
why the GCM test is not applied in Section 3.2 for modality selection) or non-
tabular (this was circumvented by using the relatively low-dimensional tabular
embedding of the chest x-ray images in Section 3.3). The PCM test, in
contrast, does not rely on this regression and is thus directly applicable in
cases where $X$ and $Z$ are high-dimensional or non-tabular modalities. The
GCM has further been adapted to settings with functional outcomes (Lundborg et
al., 2022b), continuous time stochastic processes Christgau et al. (2023),
censored outcomes (Kook et al., 2023), and extended to powerful weighted
(Scheidegger et al., 2022) and kernel-based (Fernández and Rivera, 2022)
versions. These are all COMETs proposed in the literature and we leave their
applicability in biomedical contexts as a topic for future work.

In the applications presented in this paper, random forest and LASSO
regressions were used. Random forests are computationally fast and require
little hyperparameter tuning to obtain well-performing regression estimates.
However, for very high-dimensional applications in which the number of
features exceed the number of observations, the LASSO is a fast and
computationally stable alternative.

Overall, we believe that COMETs provide a useful tool for bioinformaticians to
assess significance in applications with high-dimensional and potentially non-
tabular omics and biomedical data. The increasing familiarity of data analysts
with supervised learning methods, on which COMETs rely, help safeguard the
validity of the statistical inference. Further, the algorithm-agnostic nature
of the procedures make COMETs easily adaptable to future developments in
predictive modeling.

## Acknowledgments

We thank Niklas Pfister and David Rügamer for helpful discussions. We thank
Klemens Fröhlich and Shimeng Huang for helfpful comments on the manuscript. LK
was supported by the Swiss National Science Foundation (grant no. 214457). ARL
was supported by a research grant (0069071) from Novo Nordisk Fonden.

## References

  * Ahmed et al. [2021] K. T. Ahmed, J. Sun, S. Cheng, J. Yong, and W. Zhang.  Multi-omics data integration by generative adversarial network.  _Bioinformatics_ , 38(1):179–186, 2021.  doi: 10.1093/bioinformatics/btab608.
  * Barretina et al. [2012] J. Barretina, G. Caponigro, N. Stransky, K. Venkatesan, A. A. Margolin, S. Kim, C. J. Wilson, J. Lehár, G. V. Kryukov, D. Sonkin, et al.  The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity.  _Nature_ , 483(7391):603–607, 2012.  doi: 10.1038/nature11003.
  * Bellot and van der Schaar [2019] A. Bellot and M. van der Schaar.  Conditional independence testing using generative adversarial networks.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper_files/paper/2019/file/dc87c13749315c7217cdc4ac692e704c-Paper.pdf.
  * Berk et al. [2013] R. Berk, L. Brown, A. Buja, K. Zhang, and L. Zhao.  Valid post-selection inference.  _The Annals of Statistics_ , pages 802–837, 2013.  doi: 10.1214/12-AOS1077.
  * Berrett et al. [2019] T. B. Berrett, Y. Wang, R. F. Barber, and R. J. Samworth.  The Conditional Permutation Test for Independence While Controlling for Confounders.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(1):175–197, 10 2019.  doi: 10.1111/rssb.12340.
  * Candès et al. [2018] E. Candès, Y. Fan, L. Janson, and J. Lv.  Panning for Gold: ‘Model-X’ Knockoffs for High Dimensional Controlled Variable Selection.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 80(3):551–577, 2018.  doi: 10.1111/rssb.12265.
  * Chaudhary et al. [2018] K. Chaudhary, O. B. Poirion, L. Lu, and L. X. Garmire.  Deep learning–based multi-omics integration robustly predicts survival in liver cancer.  _Clinical Cancer Research_ , 24(6):1248–1259, 2018.  doi: 10.1158/1078-0432.CCR-17-0853.
  * Cheerla and Gevaert [2019] A. Cheerla and O. Gevaert.  Deep learning with multimodal representation for pancancer prognosis prediction.  _Bioinformatics_ , 35(14):i446–i454, 2019.  doi: 10.1093/bioinformatics/btz342.
  * Christgau et al. [2023] A. M. Christgau, L. Petersen, and N. R. Hansen.  Nonparametric conditional local independence testing.  _The Annals of Statistics_ , 51(5):2116–2144, 2023.  doi: 10.1214/23-AOS2323.
  * Fernández and Rivera [2022] T. Fernández and N. Rivera.  A general framework for the analysis of kernel-based tests.  _arXiv preprint_ , 2022.  doi: 10.48550/arXiv.2209.00124.
  * Glocker et al. [2023] B. Glocker, C. Jones, M. Roschewitz, and S. Winzeck.  Risk of bias in chest radiography deep learning foundation models.  _Radiology: Artificial Intelligence_ , 5(6):e230060, 2023.  doi: 10.1148/ryai.230060.
  * Greenland [2019] S. Greenland.  Valid p-values behave exactly as they should: Some misleading criticisms of p-values and their resolution with s-values.  _The American Statistician_ , 73(sup1):106–114, 2019.  doi: 10.1080/00031305.2018.1529625.
  * Guo and Shah [2023] F. R. Guo and R. D. Shah.  Rank-transformed subsampling: inference for multiple data splitting and exchangeable p-values.  _arXiv preprint_ , 2023.  doi: 10.48550/arXiv.2301.02739.
  * Hastie et al. [2009] T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman.  _The elements of statistical learning: data mining, inference, and prediction_ , volume 2.  Springer, 2009.
  * Johnson et al. [2019] A. E. Johnson, T. J. Pollard, N. R. Greenbaum, M. P. Lungren, C.-y. Deng, Y. Peng, Z. Lu, R. G. Mark, S. J. Berkowitz, and S. Horng.  Mimic-cxr-jpg, a large publicly available database of labeled chest radiographs.  _arXiv preprint_ , 2019.  doi: 10.48550/arXiv.1901.07042.
  * Kim et al. [2022] I. Kim, M. Neykov, S. Balakrishnan, and L. Wasserman.  Local permutation tests for conditional independence.  _The Annals of Statistics_ , 50(6):3388–3414, 2022.  doi: 10.1214/22-AOS2233.
  * Kook et al. [2023] L. Kook, S. Saengkyongam, A. R. Lundborg, T. Hothorn, and J. Peters.  Model-Based Causal Feature Selection for General Response Types.  _arXiv preprint_ , 2023.  doi: 10.48550/arXiv.2309.12833.
  * LeCun et al. [2015] Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  _Nature_ , 521(7553):436–444, 2015.  doi: 10.1038/nature14539.
  * Lundborg et al. [2022a] A. R. Lundborg, I. Kim, R. D. Shah, and R. J. Samworth.  The projected covariance measure for assumption-lean variable significance testing.  _arXiv preprint_ , 2022a.  doi: 10.48550/arXiv.2211.02039.
  * Lundborg et al. [2022b] A. R. Lundborg, R. D. Shah, and J. Peters.  Conditional independence testing in hilbert spaces with applications to functional data analysis.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 84(5):1821–1850, 2022b.  doi: 10.1111/rssb.12544.
  * Poirion et al. [2021] O. B. Poirion, Z. Jing, K. Chaudhary, S. Huang, and L. X. Garmire.  Deepprog: an ensemble of deep-learning and machine-learning models for prognosis prediction using multi-omics data.  _Genome medicine_ , 13(1):1–15, 2021.  doi: 10.1186/s13073-021-00930-x.
  * R Core Team [2021] R Core Team.  _R : A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2021.  URL https://www.R-project.org/.
  * Scheidegger et al. [2022] C. Scheidegger, J. Hörrmann, and P. Bühlmann.  The weighted generalised covariance measure.  _The Journal of Machine Learning Research_ , 23(1):12517–12584, 2022.  URL http://jmlr.org/papers/v23/21-1328.html.
  * Sellergren et al. [2022] A. B. Sellergren, C. Chen, Z. Nabulsi, Y. Li, A. Maschinot, A. Sarna, J. Huang, C. Lau, S. R. Kalidindi, M. Etemadi, et al.  Simplified transfer learning for chest radiography models using less data.  _Radiology_ , 305(2):454–465, 2022.
  * Shah and Bühlmann [2023] R. D. Shah and P. Bühlmann.  Double-estimation-friendly inference for high-dimensional misspecified models.  _Statistical Science_ , 38(1):68–91, 2023.  doi: 10.1214/22-STS850.
  * Shah and Peters [2020] R. D. Shah and J. Peters.  The Hardness of Conditional Independence Testing and the Generalised Covariance Measure.  _The Annals of Statistics_ , 48(3):1514–1538, 2020.  doi: 10.1214/19-aos1857.
  * Shi et al. [2021] C. Shi, T. Xu, W. Bergsma, and L. Li.  Double generative adversarial networks for conditional independence testing.  _The Journal of Machine Learning Research_ , 22(1):13029–13060, 2021.  URL http://jmlr.org/papers/v22/21-0294.html.
  * Smucler and Rotnitzky [2022] E. Smucler and A. Rotnitzky.  A note on efficient minimum cost adjustment sets in causal graphical models.  _Journal of Causal Inference_ , 10(1):174–189, 2022.  doi: doi:10.1515/jci-2022-0015.
  * Stahlschmidt et al. [2022] S. R. Stahlschmidt, B. Ulfenborg, and J. Synnergren.  Multimodal deep learning for biomedical data fusion: a review.  _Briefings in Bioinformatics_ , 23(2):bbab569, 2022.  doi: 10.1093/bib/bbab569.
  * Strobl et al. [2019] E. V. Strobl, K. Zhang, and S. Visweswaran.  Approximate kernel-based conditional independence tests for fast non-parametric causal discovery.  _Journal of Causal Inference_ , 7(1):20180017, 2019.  doi: 10.1515/jci-2018-0017.
  * Tay et al. [2023] J. K. Tay, B. Narasimhan, and T. Hastie.  Elastic net regularization paths for all generalized linear models.  _Journal of Statistical Software_ , 106(1):1–31, 2023.  doi: 10.18637/jss.v106.i01.
  * Williamson et al. [2021] B. D. Williamson, P. B. Gilbert, M. Carone, and N. Simon.  Nonparametric variable importance assessment using machine learning techniques.  _Biometrics_ , 77(1):9–22, 2021.  doi: 10.1111/biom.13392.
  * Williamson et al. [2023] B. D. Williamson, P. B. Gilbert, N. R. Simon, and M. Carone.  A general framework for inference on algorithm-agnostic variable importance.  _Journal of the American Statistical Association_ , 118(543):1645–1658, 2023.  doi: 10.1080/01621459.2021.2003200.
  * Wright and Ziegler [2017] M. N. Wright and A. Ziegler.  ranger: A fast implementation of random forests for high dimensional data in C++ and R.  _Journal of Statistical Software_ , 77(1):1–17, 2017.  doi: 10.18637/jss.v077.i01.
  * Zhang et al. [2012] K. Zhang, J. Peters, D. Janzing, and B. Schölkopf.  Kernel-based conditional independence test and application in causal discovery.  _arXiv preprint_ , 2012.  doi: 10.48550/arXiv.1202.3775.


# Dust-acoustic rogue waves in an electron-positron-ion-dust plasma medium

M. H. Rahman∗,1, N. A. Chowdhury†,1,2, A. Mannan‡,1,3, and A. A. Mamun§,1,4
1Department of Physics, Jahangirnagar University, Savar, Dhaka-1342,
Bangladesh
2Plasma Physics Division, Atomic Energy Centre, Dhaka-1000, Bangladesh
3Institut für Mathematik, Martin Luther Universität Halle-Wittenberg, Halle,
Germany
4Wazed Miah Science Research Center, Jahangirnagar University, Savar,
Dhaka-1342, Bangladesh
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A precise theoretical investigation has been made on dust-acoustic (DA) waves
(DAWs) in a four components dusty plasma medium having inertial warm adiabatic
dust grains and inertialess $q$-distributed electrons as well as isothermal
ions and positrons. The nonlinear and dispersive parameters of the nonlinear
Schrödinger equation (NLSE), which develops by using reductive perturbations
technique, have been used to recognize the stable and unstable parametric
regions of the DAWs as well as associated DA rogue waves (DARWs) in the
unstable parametric regime of the DAWs. The effects of the light positrons and
massive dust grains in determining the amplitude and width of the DARWs
associated with DAWs are examined. The findings of our present investigation
may be useful for understanding different nonlinear electrostatic phenomena in
both space and laboratory plasmas.

###### keywords:

Reductive perturbation method , Modulational instability , Dust-acoustic waves
, Rogue waves.

††journal: “Communications in Theoretical Physics”

##  1 Introduction

The electron-positron-ion-dust (EPID) plasma is a fully ionized gas comprises
of electrons and positrons of equal masses but opposite polarity as well as
ions and massive dust grains [1, 2, 3, 4]. The simultaneous co-existence of
the light positrons and heavy dust grains as well as electrons and ions in
contrast to the usual plasma containing of electrons and ions in space and
laboratory EPID plasma medium (EPIDPM) has rigourously changed the basic
features of the propagation of nonlinear electrostatic waves, viz., Dust-
acoustic (DA) waves (DAWs) [3, 4, 5, 6], DA solitary waves (DA-SWs) [3, 4], DA
shock waves (DA-SHWs) [7], DA rogue waves (DA-RWs) [8, 9, 10, 11], dust-ion-
acoustic waves (DI-AWs) [12, 13] in an EPIDPM, and has received much
attraction to solve the profound mystery, and has identified in space plasmas,
viz., active galactic nuclei [3, 4], pulsar magnetospheres [3, 4],
interstellar clouds [3, 4], supernova environments [3, 4], early universe
[12], the inner regions of the accretion disks surrounding the black hole as
well as in laboratory experiments of cluster explosions by intense laser beams
[14].

A new velocity distribution recognized as “Tsallis statistics” has attracted
much interest among the plasma physicists and is believed to be a useful
generalization of the conventional Boltzmann-Gibbs statistics which is not
valid for describing the long-range interaction in system such as complex
plasma system, and also has considered as an inevitable tool for analysing and
predicting the statistical features of long-range interaction in complex
system. This Tsallis distribution was first proposed by Tsallis [15] and was
first acknowledged by Renyi [16]. The degree of non-extensivity of the long-
range interaction in complex system is denoted by the entropic index $q$, and
when $(q\rightarrow 1)$, Tsallis distribution reduces to well known Maxwell-
Boltzmann velocity distribution [5, 6, 7]. Emamuddin et al. [5] examined DA
Gardner solitons in a four components dusty plasma medium (DPM) having non-
extensive plasma species. Eslami et al. [6] analyzed the properties of DA-SWs
in presence of non-extensive plasma species, and found that the DA-SWs exhibit
compression and rarefaction according to the negative and positive $q$,
respectively. Roy et al. [7] demonstrated DA-SHWs in a three components DPM
featuring non-extensive electrons, and observed that DA-SHWs potential
decreases with the increase in the value of $q$ for positive limit.

The signature of the positrons in EPIDPM has encouraged many authors to
examine the propagation of nonlinear electrostatic pulses in EPIDPM because
the effect of positrons in EPIDPM with substantial amount can not be ignored.
Banerjee and Maitra [12] considered four components EPIDPM and studied
electrostatic potential structures in presence of the massive dust grains and
light positrons and observed that the hight of the potential structures
increases with increasing value of massive dust number density but decreases
with increasing light positron number density. Paul and Bandyopadhyay [13]
demonstrated DI-AWs in an EPIDPM and highlighted the existence of both
polarities solitary waves as well as positive and negative potential double
layers but only positive super-solitons in an EPIDPM having positrons.

The nonlinear Schrödinger equation (NLSE) [8, 9, 10, 11, 17, 18, 19] has
considered one of the elegant equations which can art the picture of the
modulational instability (MI) and the mechanism of rogue waves (RWs) in any
complex medium, and also has manifested as a central archetype to investigate
the nonlinear property of complex plasma medium. Bains et al. [9] studied the
propagation of DAWs in a three components DPM having non-extensive plasma
species and found that the critical wave number ($k_{c}$) increases with
increasing the non-extensivity of the plasma species. Moslem et al. [10]
analyzed the DAWs by considering $q$-distributed plasma species and observed
that non-extensive plasma parameter $q$ plays a significant role in maximizing
or minimizing the energy of the DA-RWs. Rahman et al. [11] examined DA-RWs in
a multi-component DPM and observed that the ion temperature enhances the
height and thickness of the DA-RWs while the electron temperature reduces the
height and thickness of the DA-RWs.

Recently, Esfandyari-Kalejahi et al. [3] investigated the nonlinear
propagation of DA-SWs in an EPIDPM and observed that the amplitude of DA-SWs
increases with the increase in the value of negatively charged dust charge
state. Jehan et al. [4] studied DA-SWs in a four components DPM by considering
inertialess iso-thermal electrons, positrons, and ions as well as inertial
massive dust grains, and found that the amplitude of dip solitons increases
with increasing the value of the positron concentration. However, to the best
knowledge of the authors, there is no theoretical investigation has been made
on DAWs in EPIDPM considering warm adiabatic dust grains. Therefore, in this
paper, we are interested to study the MI of DAWs and the formation of DA-RWs
in an EPIDPM as well as the effects of $q$-distributed electrons, and
isothermal ions and positrons on the DA-RWs in a four components EPIDPM.

The rest of the paper is organized in the following fashion: The governing
equations of our plasma model are presented in Sec. 2. The MI and RWs are
given in Sec. 3. Finally, a brief conclusion is provided in Sec. 4.

##  2 Governing equations

We consider the propagation of DAWs in an unmagnetized collisionless EPIDPM
consisting of inertial warm adiabatic negatively charged massive dust grains
(mass $m_{d}$; charge $q_{d}=-Z_{d}e$), and inertialess $q$-distributed
electrons (mass $m_{e}$; charge $-e$) as well as iso-thermal ions (mass
$m_{i}$; charge $q_{i}=+Ze$) and positrons (mass $m_{p}$; charge $+e$), where
$Z_{d}$ ($Z_{i}$) is the number of electron (proton) residing on a negatively
(positively) charged massive dust grains (ions). Overall, the charge
neutrality condition for our plasma model can be written as
$n_{e0}+Z_{d}n_{d0}=n_{p0}+Z_{i}n_{i0}$. Now, the basic set of normalized
equations can be written in the form

$\displaystyle\frac{\partial n_{d}}{\partial
t}+\frac{\partial(n_{d}u_{d})}{\partial x}=0,$ (1)
$\displaystyle\frac{\partial u_{d}}{\partial t}+u_{d}\frac{\partial
u_{d}}{\partial x}+\delta n_{d}\frac{\partial n_{d}}{\partial
x}=\frac{\partial\varphi}{\partial x},$ (2)
$\displaystyle\frac{\partial^{2}\varphi}{\partial
x^{2}}=(\mu_{p}+\mu_{i}-1)n_{e}-\mu_{p}n_{p}-\mu_{i}n_{i}+n_{d},$ (3)

where $n_{d}$ is the adiabatic dust grains number density normalized by its
equilibrium value $n_{d0}$; $u_{d}$ is the dust fluid speed normalized by the
DA wave speed $C_{d}=(Z_{d}K_{B}T_{i}/m_{d})^{1/2}$ (with $T_{i}$ being the
ion temperature, $m_{d}$ being the dust grain mass, and $K_{B}$ being the
Boltzmann constant); $\phi$ is the electrostatic wave potential normalized by
$K_{B}T_{i}/e$ (with $e$ being the magnitude of single electron charge); the
time and space variables are normalized by $\omega_{pd}^{-1}=(m_{d}/4\pi
Z_{d}^{2}e^{2}n_{d0})^{1/2}$ and $\lambda_{Dd}=(K_{B}T_{i}/4\pi
Z_{d}n_{d0}e^{2})^{1/2}$, respectively; $P_{d}=P_{d0}(N_{d}/n_{d0})^{\gamma}$
[with $P_{d0}$ being the equilibrium adiabatic pressure of the dust, and
$\gamma=(N+2)/N$, where $N$ is the degree of freedom and for one-dimensional
case, $N=1$ then $\gamma=3$]; $P_{d0}=n_{d0}K_{B}T_{d}$ (with $T_{d}$ being
the temperatures of the adiabatic dust grains); and other plasma parameters
are considered as $\delta=3T_{d}/Z_{d}T_{i}$, $\mu_{p}=n_{p0}/Z_{d}n_{d0}$,
and $\mu_{i}=Z_{i}n_{i0}/Z_{d}n_{d0}$. Now, the expression for electron number
density obeying $q$-distribution is given by [15, 8]

$\displaystyle n_{e}=\left[1+(q-1)\sigma\varphi\right]^{\frac{q+1}{2(q-1)}},$
(4)

where $\sigma=T_{i}/T_{e}$ and $T_{e}>T_{i}$. The expression for the number
density of ions and positrons obeying iso-thermal distribution are,
respectively, given by

$\displaystyle n_{p}=\mbox{exp}(-\alpha\varphi),$ (5) $\displaystyle
n_{i}=\mbox{exp}(-\varphi),$ (6)

where $\alpha=T_{i}/T_{p}$ and $T_{p}>T_{i}$. Now, by substituting Eqs.
(4)-(6) in Eq. (3), and expanding up to third order in $\varphi$, we get

$\displaystyle\frac{\partial^{2}\varphi}{\partial
x^{2}}+1=n_{d}+\Upsilon_{1}\varphi+\Upsilon_{2}\varphi^{2}+\Upsilon_{3}\varphi^{3}+\cdot\cdot\cdot,$
(7)

where

$\displaystyle\Upsilon_{1}=[\sigma(q+1)(\mu_{p}+\mu_{i}-1)+2\mu_{p}\alpha+2\mu_{i}]/2,$
$\displaystyle\Upsilon_{2}=[\sigma^{2}(q+1)(q-3)(1-\mu_{p}-\mu_{i})-4\mu_{p}\alpha^{2}-4\mu_{i}]/8,$
$\displaystyle\Upsilon_{3}=[\sigma^{3}(q+1)(q-3)(3q-5)(\mu_{p}+\mu_{i}-1)+F1]/48,\
$

where $F1=8\mu_{p}\alpha^{3}+8\mu_{i}$. To study the MI of the DAWs, we want
to derive the NLSE by employing the reductive perturbation method and for that
case, first we can write the stretched coordinates as
$\xi={\varepsilon}(x-v_{g}t)$ and $\tau={\varepsilon}^{2}t$, where
$\varepsilon$ is a smallness parameter and $v_{g}$ is the group velocity of
DAWs to be determined later. Now, we can write the dependent variables as [20,
21, 22]

$\displaystyle
n_{d}=1+\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}n_{dl}^{(m)}(\xi,\tau)~{}\mbox{exp}[il(kx-\omega
t)],$ (8) $\displaystyle
u_{d}=\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}u_{dl}^{(m)}(\xi,\tau)~{}\mbox{exp}[il(kx-\omega
t)],$ (9)
$\displaystyle\varphi=\sum_{m=1}^{\infty}\epsilon^{m}\sum_{l=-\infty}^{\infty}\varphi_{l}^{(m)}(\xi,\tau)~{}\mbox{exp}[il(kx-\omega
t)],$ (10)

where $k$ and $\omega$ is real variables presenting the carrier wave number
and frequency, respectively. The derivative operators can be written as

$\displaystyle\frac{\partial}{\partial t}\rightarrow\frac{\partial}{\partial
t}-\epsilon
v_{g}\frac{\partial}{\partial\xi}+\epsilon^{2}\frac{\partial}{\partial\tau},$
(11) $\displaystyle\frac{\partial}{\partial
x}\rightarrow\frac{\partial}{\partial
x}+\epsilon\frac{\partial}{\partial\xi}.$ (12)

Now, by substituting Eqs. (8)-(12) into Eqs. (1), (2), and (7), and taking the
terms containing $\varepsilon$, the first-order ($m=1$ with $l=1$) reduced
equation can be written as

$\displaystyle
n_{d1}^{(1)}=\frac{k^{2}}{\Omega}\varphi_{1}^{(1)},~{}~{}~{}~{}~{}~{}~{}u_{d1}^{(1)}=\frac{k\omega}{\Omega}\varphi_{1}^{(1)},$
(13)

where $\Omega=\delta k^{2}-\omega^{2}$. We thus obtain the dispersion relation
for DAWs

$\displaystyle\omega^{2}=\frac{k^{2}}{\Upsilon_{1}+k^{2}}+\delta k^{2},$ (14)

The second-order ($m=2$ with $l=1$) equations are given by

$\displaystyle
n_{d1}^{(2)}=\frac{k^{2}}{\Omega}\varphi_{1}^{(2)}-\frac{2ik\omega(v_{g}k-\omega)}{\Omega^{2}}\frac{\partial\varphi_{1}^{(1)}}{\partial\xi},$
(15) $\displaystyle
u_{d1}^{(2)}=\frac{k\omega}{\Omega}\varphi_{1}^{(2)}-\frac{i(v_{g}k-\omega)(\omega^{2}+\delta
k^{2})}{\Omega^{2}}\frac{\partial\varphi_{1}^{(1)}}{\partial\xi},$ (16)

with the compatibility condition

$\displaystyle v_{g}=\frac{\partial\omega}{\partial
k}=\frac{\omega^{2}-\Omega^{2}}{k\omega}.$ (17)

The coefficients of $\epsilon$ for $m=2$ and $l=2$ provide the second order
harmonic amplitudes which are found to be proportional to
$|\varphi_{1}^{(1)}|^{2}$

$\displaystyle n_{d2}^{(2)}=B_{1}|\varphi_{1}^{(1)}|^{2},$ (18) $\displaystyle
u_{d2}^{(2)}=B_{2}|\varphi_{1}^{(1)}|^{2},$ (19)
$\displaystyle\varphi_{2}^{(2)}=B_{3}|\varphi_{1}^{(1)}|^{2},$ (20)

where

$\displaystyle B_{1}=\frac{2B_{3}k^{2}\Omega^{2}-3\omega^{2}k^{4}-\delta
k^{6}}{2\Omega^{3}},$ $\displaystyle B_{2}=\frac{2B_{1}\delta
k\Omega^{2}+\delta
k^{5}+\omega^{2}k^{3}-2B_{3}k\Omega^{2}}{2\omega\Omega^{2}},$ $\displaystyle
B_{3}=\frac{3\omega^{2}k^{4}+\delta
k^{6}-2\Upsilon_{2}\Omega^{3}}{6k^{2}\Omega^{3}}.\ $

Now, we consider the expression for ($m=3$ with $l=0$) and ($m=2$ with $l=0$),
which leads the zeroth harmonic modes. Thus, we obtain

$\displaystyle n_{d0}^{(2)}=B_{4}|\varphi_{1}^{(1)}|^{2},$ (21) $\displaystyle
u_{d0}^{(2)}=B_{5}|\varphi_{1}^{(1)}|^{2},$ (22)
$\displaystyle\varphi_{0}^{(2)}=B_{6}|\varphi_{1}^{(1)}|^{2},$ (23)

where

$\displaystyle B_{4}=\frac{B_{6}\Omega^{2}-2v_{g}\omega k^{3}-\delta
k^{4}-k^{2}\omega^{2}}{\Omega^{2}(\delta-v_{g}^{2})},$ $\displaystyle
B_{5}=\frac{B_{4}\delta\Omega^{2}-B_{6}\Omega^{2}+\delta
k^{4}+k^{2}\omega^{2}}{v_{g}\Omega^{2}},$ $\displaystyle
B_{6}=\frac{2v_{g}\omega k^{3}+\delta
k^{4}+k^{2}\omega^{2}-2\Upsilon_{2}\Omega^{2}(\delta-
v_{g}^{2})}{\Omega^{2}(1+\Upsilon_{1}\delta-\Upsilon_{1}v^{2}_{g})}.\ $

Finally, the third harmonic modes ($m=3$) and ($l=1$), with the help of
(13)$-$(23), give a set of equations, which can be reduced to the following
NLSE:

$\displaystyle
i\frac{\partial\phi}{\partial\tau}+P\frac{\partial^{2}\phi}{\partial\xi^{2}}+Q|\phi|^{2}\phi=0,$
(24)

where $\phi=\varphi_{1}^{(1)}$ for simplicity. In Eq. (24), $P$ is the
dispersion coefficient which can be written as

$\displaystyle
P=\frac{(v_{g}k-\omega)(\omega^{3}-3v_{g}k\omega^{2}+3\delta\omega
k^{2}-v_{g}\delta k^{3})-\Omega^{3}}{2\omega\Omega k^{2}},\ $

and $Q$ is the nonlinear coefficient which can be written as

$\displaystyle
Q=\frac{3\Upsilon_{3}\Omega^{2}+2\Upsilon_{2}\Omega^{2}(B_{3}+B_{6})-2\omega
k^{3}(B_{2}+B_{5})-F2}{2\omega k^{2}},\ $

where $F2=(\delta k^{4}+k^{2}\omega^{2})(B_{1}+B_{4})$. The space and time
evolution of the DAWs in EPIDPM are directly governed by the dispersion ($P$)
and nonlinear ($Q$) coefficients, and indirectly governed by different plasma
parameters such as as $\delta$, $\mu_{p}$, $\mu_{i}$, $\sigma$, $\alpha$, $q$,
and $k$. Thus, these plasma parameters significantly affect the stability
conditions of the DAWs.

Figure 1: $P/Q$ versus $k$ curve for different values of positive $q$ along
with $\alpha=0.3$, $\delta=0.005$, $\mu_{i}=1.4$, $\mu_{p}=0.3$, and
$\sigma=0.7$. Figure 2: $P/Q$ versus $k$ curve for different values of
negative $q$ along with $\alpha=0.3$, $\delta=0.005$, $\mu_{i}=1.4$,
$\mu_{p}=0.3$, and $\sigma=0.7$. Figure 3: $P/Q$ versus $k$ curve for
different values of $\mu_{i}$ along with $\alpha=0.3$, $\delta=0.005$,
$\mu_{p}=0.3$, $\sigma=0.7$ and $q=1.5$. Figure 4: $P/Q$ versus $k$ curve for
different values of $\alpha$ along with $\delta=0.005$, $\mu_{i}=1.4$,
$\mu_{p}=0.3$, $\sigma=0.7$, and $q=1.5$. Figure 5: variation of $|\phi|$ with
$\xi$ for different values of $\mu_{p}$ along with $\alpha=0.3$,
$\delta=0.005$, $\mu_{i}=1.4$, $\sigma=0.7$, and $q=1.5$. Figure 6: variation
of $|\phi|$ with $\xi$ for different values of $\sigma$ along with
$\alpha=0.3$, $\delta=0.005$, $\mu_{i}=1.4$, $\mu_{p}=0.3$, and $q=1.5$.

##  3 Modulational instability and rogue waves

The stability of DAWs in an EPIDPM is governed by the sign of dispersion
coefficient ($P$) and nonlinear coefficient ($Q$) which are function of
various plasma parameters such as $\alpha$, $\delta$, $\mu_{i}$, $\mu_{p}$,
$\sigma$, and $q$. These plasma parameters significantly control the stability
conditions of the DAWs. When $P$ and $Q$ are same sign ($P/Q>0$), the
evolution of the DAWs amplitude is modulationally unstable and in this region
electrostatic bright envelope solitons as well as RWs exist [8, 23]. On the
other hand, when $P$ and $Q$ are opposite sign ($P/Q<0$), the DAWs are
modulationally stable in presence of external perturbations and in this region
electrostatic dark envelope solitons exist [8, 23]. The plot of $P/Q$ against
$k$ yields stable and unstable regions of the DAWs. The point, at which
transition of $P/Q$ curve intersects with $k$-axis, is known as the threshold
or critical wave number $k$ ($=k_{c}$).

We have depicted $P/Q$ versus $k$ graph for different values of positive $q$
and negative $q$ in Figs. 1 and 2, respectively, and it is obvious from these
figures that (a) the DAWs remain stable for small $k$ ($k<k_{c}$) and the MI
sets in for values of $k$ ($k>k_{c}$); (b) when $q=1.1$, $1.5$, and $1.9$ then
the corresponding $k_{c}$ value is $k_{c}\equiv 1.88$ (dotted blue curve),
$k_{c}\equiv 1.92$ (dashed green curve), and $k_{c}\equiv 1.95$ (solid red
curve); (c) $k_{c}$ increases with increasing value of positive $q$. On the
other hand, in the case of negative $q$: (d) when $q=-0.9$, $-0.5$, and $-0.1$
then the corresponding $k_{c}$ value is $k_{c}\equiv 1.57$ (dotted blue
curve), $k_{c}\equiv 1.67$ (dashed green curve), and $k_{c}\equiv 1.72$ (solid
red curve); (e) the negative $q$ also increases the value of critical wave
number. Hence one can say that the stability of the DAWs is independent on the
sign of $q$ but dependent on the magnitude of $q$.

We have investigated the effects of ion and dust number density as well as
their charge state to organize the stable and unstable parametric regions of
DAWs by depicting $P/Q$ versus $k$ graph for different values of $\mu_{i}$ in
Fig. 3. It can be seen from this figure that (a) the $k_{c}$ increases with
$\mu_{i}$; the stable (unstable) parametric region increases (decreases) with
increasing value of ion (dust) number density for the fixed values of $Z_{d}$
and $Z_{i}$; (b) the $k_{c}$ increases with the increase in the value of
positive ion charge state ($Z_{i}$) while decreases with the increase in the
value of negative dust charge state ($Z_{d}$) for a constant value of $n_{i0}$
and $n_{d0}$ (via $\mu_{i}=Z_{i}n_{i0}/Z_{d}n_{d0}$). So, the charge and the
number density of the positive ion and negative dust grains play an opposite
role in recognizing the stable and unstable parametric regions of the DAWs.

Figure 4 represents the stability nature of the DAWs according to the ion to
positron temperature (via $\alpha=T_{i}/T_{p}$). The unstable window increases
with an increase in ion temperature while decreases with an increase in the
positron temperature.

The governing equation for the highly energetic DA-RWs in the unstable region
($P/Q>0$) can be written as [24, 25]

$\displaystyle\phi(\xi,\tau)=\sqrt{\frac{2P}{Q}}\left[\frac{4(1+4iP\tau)}{1+16P^{2}\tau^{2}+4\xi^{2}}-1\right]\mbox{exp}(2iP\tau).$
(25)

We have also numerically analyzed Eq. (25) in Fig. 5 to illustrate the
influence of number density ratio of positron to dust (via
$\mu_{p}=n_{p0}/Z_{d}n_{d0}$) on the formation of DA-RWs associated with
unstable parametric region (i.e., $P/Q>0$), and it is clear from this figure
that (a) the nonlinearity as well as the amplitude and width of the DA-RWs
increases with increasing positron number density $n_{p0}$ for a constant
value of $Z_{d}$ and $n_{d0}$; (b) On the other hand, the nonlinearity as well
as the amplitude and width of the DA-RWs decreases with increasing negative
dust number density $n_{d0}$ as well as negative dust charge state $Z_{d}$ for
a constant value of $n_{p0}$.

The effect of ion temperature ($T_{i}$) as well as electron temperature
($T_{e}$) on the height and thickness of DA-RWs can be observed from Fig. 6
and it can be manifested from this figure that the height and thickness of the
DA-RWs increase as we increase the ion temperature while decrease with an
increase in the value of electron temperature; (b) physically, the
nonlinearity of the plasma medium enhances with ion temperature by developing
a gigantic DA-RWs associated with DAWs in unstable parametric region while the
nonlinearity of the plasma medium reduces with electron temperature by
developing a smaller DA-RWs associated with DAWs in unstable parametric
region. This result is a nice agreement with Rahman et al. [8] work.

##  4 Conclusion

In this study, we have performed a nonlinear analysis of DAWs in an EPIDPM
having inertial massive dust grains and inertialess $q$-distributed electrons
as well as iso-thermal positrons and ions. The nonlinear properties of the
plasma medium as well as the mechanism to generate gigantic DA-RWs associated
with DAWs are governed by the standard NLSE. It has been seen that iso-thermal
ions and positrons as well as non-extensive electrons have essential influence
on the MI of DAWs. The dependency of the height and thickness of the DA-RWs on
the various plasma parameters, namely, charge state of the negatively charged
dust, number density of the positron and dust grains as well as the
temperature of the ion and electron is also examined. To conclude, the results
of our present investigation might be useful to understand the nonlinear
phenomena (viz., MI of DAWs and DA-RWs) in space dusty plasmas, viz., active
galactic nuclei [3, 4], pulsar magnetospheres [3, 4], interstellar clouds [3,
4], supernova environments [3, 4], our early universe [12], the inner regions
of the accretion disks surrounding the black hole as well as in laboratory
experiments of cluster explosions by intense laser beams [14].

## References

  * [1] R. K. Shikha, et al., Eur. Phys. J. D 73, 177 (2019).
  * [2] S. Jahan, et al., Commun. Theor. Phys. 71, 327 (2019).
  * [3] A. Esfandyari-Kalejahi, et al., Phys. Plasmas 19, 082308 (2012).
  * [4] N. Jehan, W Masood, and A. M. Mirza, Phys. Scr. 80, 035506 (2009).
  * [5] M. Emamuddin, et al., Phys. Plasmas 20, 043705 (2013).
  * [6] P Eslami, et al., Phys. Plasmas 18, 102303 (2011).
  * [7] K. Roy, et al., Astrophys. Space Sci. 350, 599 (2014).
  * [8] M. H. Rahman, et al., Chinese J. Phys. 56, 645 (2018).
  * [9] A. S. Bains, et al., Astrophys. Space Sci. 343, 621 (2013).
  * [10] W. M. Moslem, et al., Phys. Rev. E 84, 066402 (2011).
  * [11] M. H. Rahman, et al., Phys. Plasmas 25, 102118 (2018).
  * [12] G. Banerjee and S. Maitra, Phys. Plasmas 23, 123701 (2016).
  * [13] A. Paul and A. Bandyopadhyay, Astrophys. Space Sci. 361, 172 (2016).
  * [14] C. Higdon, et al., Astrophys. J. 698, 350 (2009).
  * [15] C. Tsallis, J. Stat. Phys. 52, 479 (1988).
  * [16] A. Renyi, Acta Math. Acad. Sci. Hung. 6, 285 (1995).
  * [17] N. A. Chowdhury, et al., Chaos 27, 093105 (2017).
  * [18] N. A. Chowdhury, et al., Plasma Phys. Rep. 45, 459 (2019).
  * [19] N. Ahmed, et al., Chaos 28, 123107 (2018).
  * [20] N. A. Chowdhury, et al., Vacuum 147, 31 (2018).
  * [21] M. Hassan, et al., Commun. Theor. Phys. 71, 1017 (2019).
  * [22] N. A. Chowdhury, et al., Contrib. Plasma Phys. 58, 870 (2018).
  * [23] R. Fedele, Phys. Scr. 65, 502 (2002).
  * [24] N. Akhmediev, et al., Phys. Rev. E 80, 026601 (2009).
  * [25] A. Anikiewicz, et al., Phys. Lett. A 373, 3997 (2009).