Dynamical analysis of a predator-prey system with prey vigilance and hunting cooperation in predators

1.   Introduction

In practical applications, due to the interference of various factors, collected data is often incomplete. Missing data is common in public opinion polls, medical research, experimental science, and other application fields. Missing data will not only result in the reduction of effective information, the deviation of the estimation result, but also affect the statistical decision-making and distort the analysis result to some extent. One approach to deal with missing data is complete-case analysis, which deletes all incomplete data. However, Little and Rubin ^[1] pointed out that this will cause biased estimation when the occurrence of missing data is not completely at random. Yates ^[2] introduced an imputation method which is widely used to handle missing responses. The purpose of this method is to find suitable values for the missing data to impute. Then, the data of the filled values are regarded as the complete observation data, which can be analyzed by the classical method. Inverse probability weighting (IPW), which was proposed by Horvitz and Thompson ^[3], is another method to deal with missing data. The inverse of the selection probability is chosen to be the weight assigned to the fully observed data. The missing at random (MAR) assumption, in the sense of Rubin et al. ^[4], is a common assumption for statistical analysis with missing data.

In the case of missing data, the missing mechanism is usually unknown, and parameter methods and nonparametric methods are commonly used to estimate. For the parameter method, there may be a model misspecification problem. Imai and Ratkovic ^[5] proposed the covariate balanced propensity score (CBPS), which improves the parameter method. Based on the CBPS method, Guo et al. ^[6] applied the CBPS method to mean regression to obtain the estimators of the regression parameters $\beta$ and the mean $\mu$ in the case of missing data.

Expectile regression, which was proposed by Newey and Powell ^[7], can be regarded as a generalization of mean regression. Expectile regression uses the sum of asymmetric residual squares as the loss function, and since the loss function is convex and differentiable, expectile regression has computational advantages over quantile regression. Recently, people have carried out a lot of specific research on expectile regression. Sobotka et al. ^[8] established the asymptotic properties of a semi-parametric expectile regression estimator and introduced confidence intervals for expectiles. Waltrup et al. ^[9] observed that expectile regression tends to have less crossing and more robustness against heavy tailed distributions than quantile regression. Ziegel ^[10] concluded that expectile shares coherence and elicitability. Pan et al. ^[11] considered fitting a linear expectile regression model for estimating conditional expectiles based on a large quantity of data with covariates missing at random. Recently, Pan et al. ^[12] developed a weighted expectile regression approach for estimating the conditional expectile when covariates are missing at random (MAR). They only considered a single expectile, and the missing mechanism was assumed to be logistic regression. However, the missing mechanism model may be misspecified. In addition, it is known that making full use of multiple target information can improve the efficiency of parameter estimation. In summary, when the model may be misspecified, we use the idea of covariate balance to study the weighted expectile average estimation of unknown parameters based on CBPS by using multiple expected information. Our estimators can improve performance of the usual weighted expectile average estimator in terms of standard deviation (SD) and mean squared error (MSE).

The rest of this paper is organized as follows. In Section 2, we propose a CBPS-based estimator for the propensity score. In Section 3, we estimate the expected quantile weighted average of the regression parameters based on CBPS. Moreover, we establish the asymptotic normality of the weighted estimator in Section 4. In Section 5, a simulation study is carried out to assess the performance of the proposed method. The proofs of those theoretical results are deferred to the Appendix.

2.   CBPS-based estimator for the propensity score

Consider the following linear regression model:

where $Y_i$ is response, $X_i$ is covariate, $\beta$ is the $p$-dimensional vector of unknown parameters, and $\varepsilon_i$ is the random error. Assuming that the response variable $Y_i$ is missing at random, the covariate $X_i$ can be fully observed. For the $i^{th}$ individual, let $\delta_i$ denote the observing indicator, i.e., $\delta_i = 1$ if $Y_i$ is observed and $0$ otherwise. In our paper, we only consider the missing mechanism of missing at random (MAR), that is,

where $\pi_i$ is called the selection probability function or the propensity score.

The most popular choice of $\pi(X_i)$ is a logistic regression function (Peng et al. ^[13]). We make the same choice and posit a logistic regression model for $\pi(X_i)$,

and $\gamma = (\gamma_0, \gamma_1^T)^T\in\Theta$ is the unknown parameter vector with the parameter space $\Theta\subseteq {R^{q+1}}$. Here, $\gamma$ can be estimated by maximizing the log-likelihood function

Assuming that $\pi(X_i, \gamma)$ is twice continuously differentiable with respect to $\gamma$, maximizing $L(\gamma)$ implies the first-order condition

where $\pi'(X_i, \gamma) = \partial \pi(X_i, \gamma)/\partial \gamma^T.$ The maximum likelihood method is a commonly used and simple parameter estimation method. However, when the selection probability model (2.3) is assumed to be wrong, the estimator based on this method will have a large deviation. In order to make the parameter method more robust, we use the covariate balanced propensity score method proposed by Imai and Ratkovic ^[5] to estimate the unknown parameter $\gamma$, that is,

$\tilde{X_i} = f(X_i)$ is an M-dimensional vector-valued measurable function of $X_i$. For any covariate function, as long as the expectation exists, Eq (2.5) must hold. If the propensity score model is incorrectly specified, then the maximum likelihood may not be able to balance the covariates. Following Imai and Ratkovic ^[5], we can set $\tilde{X_i} = X_i$ to ensure that the first moment of each covariate is balanced even when the model is misspecified. $\pi(X_i, \gamma)$ satisfies the condition

The sample form of the covariate equilibrium condition obtained from (2.6) is

where

According to Imai and Ratkovic ^[5], if we only use the condition of the $\pi'(X_i, \gamma)$ equilibrium, i.e., (2.4), at this time, the number of equations is equal to the number of parameters. Then, the covariate equilibrium propensity score is just-identified. If we combine Eq (2.4) with the score condition given in Eq (2.7),

where

then the covariate equilibrium propensity score is over-identified because the number of moment conditions exceeds that of parameters. For over-identified CBPS, the estimation of $\gamma$ can be obtained by using the generalized moment method (GMM) (Hansen ^[14]). For a positive semidefinite symmetric weight matrix $W$, the GMM estimator $\hat{\gamma}$ can be obtained by minimizing the following objective function for $\gamma$:

The above method is also applicable to the case where the covariate balanced propensity score is just-identified.

3.   Estimator for the regression parameter

Pan et al. ^[12] introduced the weighted expectile regression estimation of a linear model in detail. According to the idea of inverse probability weighting, when the selection probability function $(\pi_1\dots, \pi_n)^T$ is known, the expectile estimator of $\beta$ under missing responses is defined as

where $\tau_k\in (0, 1)$ is expectile level, and $\Phi_{\tau_k}(v) = \left|\tau_k-I(v\le0) \right|v^2$. $b_{\tau_k}$ represents the ${\tau_k}$-expectile of the error term $\varepsilon_i$. Then, according to Zhao et al.^[15], let $K$ be the number of expectiles, and consider the equally spaced expectiles $\tau_k = \frac{k}{K+1}$, $k = 1, 2, \dots, K$. The weighted expectile average estimator of the linear model parameter $\beta$ when the missing mechanism is known is defined as

where the weight vector $(\omega_1, \dots, \omega_K)^T$ satisfies $\sum_{k = 1}^{K}\omega_k = 1$.

When the selection probability function is unknown, we use the method proposed in the second section to estimate the parameter $\gamma$ based on CBPS, so as to obtain $\pi(X_i, \hat{\gamma})$. The loss function of the $\tau_k$-expectile can be defined as

By minimizing the loss function, we can obtain the expectile estimation of the unknown parameter $\beta$,

Therefore, the weighted expectile average estimation of the linear model parameter $\beta$ when the missing mechanism is unknown under the missing responses is defined as

The weight vector $(\omega_1, \dots, \omega_K)^T$ satisfies $\sum_{k = 1}^{K}\omega_k = 1$.

4.   Asymptotic property

Let $\gamma_0$ and $\beta_0$ represent the true values of $\gamma$ and $\beta$ respectively, and $U(\gamma) = \begin{pmatrix} s(\delta, X, \gamma)\\z(\delta, X, \gamma)X\end{pmatrix}$. In addition, with reference to Pan et al. ^[12] and Guo ^[16], the following regularity conditions are required.

C1: $\gamma_0$ is the interior point of $\Theta$.

C2: $U(\gamma)$ is differentiable in the neighborhood $\triangle$ of $\gamma_0$.

C3: $E[U(\gamma_0)] = 0$, $E[{\left \| U(\gamma_0) \right\|}^2] < \infty$.

C4: $E[\sup_{\gamma\in \triangle }\left \| \nabla_\gamma U(\gamma) \right \|] < \infty$, where $\nabla_{\gamma}$ is the first-order partial derivative of the function to $\gamma$.

C5: $\Gamma = E[\nabla_\gamma U(\gamma)]$ exists.

C6: For any $i$, there exists a compact set $\mathcal{X}$, such that $X_i\in\mathcal{X}\subset\mathbb{R}^p$, and $X_i$ and $\varepsilon_i$ are independent.

C7: The regression errors $\left\{\varepsilon_i\right\}_{i = 1}^n$ are independent and identically distributed with common cumulative distribution function $F(\cdot)$, satisfying $E[\varepsilon_i^2] < \infty$.

C8: There exists $a > 0$ such that $\pi(V_i, \gamma) > a$ for any $i$.

C9: The symmetric matrix $\Sigma_1$ is positive definite.

The following theorem presents the asymptotic distribution of $\hat{\beta}_w$.

Theorem 4.1 (Asymptotic Normality of $\hat{\beta}_w$) Under the assumptions C1–C9, we have

where $\Sigma_1 = E[X_iX_i^T]$, $\Lambda = E[\lambda\lambda^T], \lambda = \mu-E[\partial \mu/\partial \gamma^T]\{E[\partial U(\gamma)/\partial \gamma^T]\}^{-1}U(\gamma)$, $\mu = \frac{\delta}{\pi(X, \gamma)}X\sum_{k = 1}^{K}\frac{\omega_k\Psi_{\tau_k}(\varepsilon-b_{0k})}{g(\tau_k)}$.

5.   Simulation

In the following, the expectile weighted average estimation based on covariate balancing propensity score proposed in this paper is analyzed by numerical simulation, and the method is compared with the usual parameter estimation method in the case of correct and wrong model assumptions. Consider the following linear model:

where $\beta_1 = 0.5$, $\beta_2 = 1$, $\beta_3 = 1$, and $(X_1, X_2, X_3)$ obeys the joint normal distribution with mean of $0$, covariance of $0.5$, and variance of $1$. The error term $\varepsilon$ obeys the standard normal distribution. In our simulation, we take $K = 10$, $\tau_{k} = k/11$ for $k = 1, 2, \dots, 10$, and consider the real choice probability model as

Under the assumption of random missing, in order to illustrate the effect when the model is misspecified, we assume that the covariates

If the model (5.2) is represented by $\pi(X^*)$, the model will be specified incorrectly. In the simulation study of the expectile regression of the unknown parameter $\beta$, we consider the following two cases: (1) Propensity score model is correctly specified. (2) Propensity score model is misspecified. Zhao ^[15] proposed the weighted composite expectile regression method for a varying-coefficient partially linear model. For a given scenario, referring to Zhao ^[15], we compare the weighted expectile average estimation based on CBPS, denoted as CBPS-WEAE, with weighted composite expectile regression, denoted as WCER, and weighted composite quantile regression, denoted as WCQR, to examine the performance of the estimator, where the weights of WCER and WCQR are estimated by the generalized linear model.

In the simulation, samples of size $n = 500,800, 1000, 1200$ are generated independently. For each scenario, we conduct $1000$ simulations and calculate the average mean squared error (MSE) for estimator of $\beta$ and the average bias (Bias) and standard deviation (SD) for estimator of $\beta_1$, $\beta_2$, and $\beta_3$. In order to examine the influence of the error distribution on the performance of the proposed method, two different distributions of the model error $\varepsilon$ are considered: standard normal distribution $N(0, 1)$ and centralized $\chi^2$ distribution with 4 degrees of freedom. The results of our simulations are presented in Tables 1 and 2.

From Tables 1 and 2 we observe that, as expected, all three estimators are unbiased. In terms of MSE, as a convenient measure of average error, we observe that when model error $\varepsilon$ follows the standard normal distribution $N(0, 1)$, CBPS-WEAE performs best among the three estimators considered, followed immediately by WCER, while WCQR performs worst. When $\varepsilon$ follows a centralized $\chi^2$ distribution with 4 degrees of freedom, CBPS-WEAE is superior to the other two methods. When sample size is large, it can be seen that the performance of the three estimators is significantly improved compared with that when the sample size is small. In general, our proposed improved estimator is effective.

6.   Conclusions

In this paper, in order to improve the estimation efficiency of weighted expectile average estimation, we estimate the selection probability function based on CBPS and propose a weighted expectile average estimator based on CBPS when the response variables are missing at random. The asymptotic normality of the proposed method is proved, and the estimation effect of the method is further illustrated by numerical simulation. The numerical simulation results show that the method is effective.

Author contributions

Qiang Zhao: Conceptualization, methodology, supervision, writing-review and editing; Zhaodi Wang: Validation, software, writing-original draft; Jingjing Wu: Funding acquisition, formal analysis, writing-original draft; Xiuli Wang: Funding acquisition, investigation, resources, writing-review and editing. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgments

The research is supported by Natural Science Foundation of Shandong Province (Grant Nos. ZR2021MA077 and ZR2021MA048).

Conflict of interest

All authors declare that there is no conflict of interest.

Appendix: Assumptions and proofs

Define the following symbols:

$\eta_i = \frac{\delta_i}{\pi(X_i, \gamma)}X_i\Psi_{\tau_k}(\varepsilon_i)$,

$\hat{\eta_i} = \frac{\delta_i}{\pi(X_i, \hat{\gamma})}X_i\Psi_{\tau_k}(\varepsilon_i)$,

$F_n = -\frac{1}{\sqrt{n}}\sum_{i = 1}^{n}\hat{\eta}_i$,

$\varepsilon_i = Y_i-X_i^T\beta_0$,

$\omega = (\omega_1, \omega_2, ...\omega_n)^T$,

$\Sigma_1 = E[X_iX_i^T]$,

$\Psi_{\tau_k} = 2\left|\tau_k-I(v\le0) \right|v$,

$u = (u_1, u_2, ..., u_n)^T$,

$G_n(u) = \sum_{i = 1}^{n}\frac{\delta_i}{\pi(X_i, \hat{\gamma})}\left[\Psi_{\tau_k}(\varepsilon_i-\frac{X_i^Tu}{\sqrt{n}})-\Psi_{\tau_k}(\varepsilon_i)\right]$.

Lemma 1. Assume that C1–C5 hold. Then, when $n\to \infty$,

where $\Gamma = E[\nabla_{\gamma}U(\gamma)]$, $\Sigma = E[U(\gamma)U^T(\gamma)]$.

The proof of Lemma 1 can refer to Theorem 2.2.1 in Guo ^[16].

Lemma 2. If the conditions C1–C4 are satisfied, then

where $\Omega = E[QQ^T], Q = \eta-E[\partial \eta/\partial \gamma^T]\{E[\partial U(\gamma)/\partial \gamma^T]\}^{-1}U(\gamma)$.

Proof. By expanding $\frac{1}{\sqrt{n}}\sum_{i = 1}^{n}\hat{\eta}_i$ at $\gamma$ and the proof process of Lemma 1, we can get

where $D_n = \left[\frac{1}{n}\sum_{i = 1}^{n}\frac{\partial \eta_i}{\partial \gamma} \right]_{\gamma^\ast}, B_n = \left[\frac{1}{n}\sum_{i = 1}^{n}\frac{\partial U_i(\gamma)}{\partial \gamma}\right]_{\gamma^\ast}$, and $\gamma^\ast$ lies between $\gamma$ and $\hat{\gamma}$.

According to the central limit theorem,

where $\Omega = E[QQ^T], Q = \eta-E[\partial \eta/\partial \gamma^T]\{E[\partial U(\gamma)/\partial \gamma^T]\}^{-1}U(\gamma)$.

Therefore, Lemma 2 is proved.

Lemma 3. If the conditions C1–C4 are satisfied, then

Proof. If the conditions C1–C4 are satisfied, it can be known from Pan et al. ^[12] that

where $g(\tau) = (1-\tau)F(0)+\tau(1-F(0))$.

Known by Hjort and Pollard ^[17], if

where $D_n(u)$ is a convex objective function with minimum point $\hat{u}_n$, $A$ is a symmetric and positive definite matrix, and $B$ is a random variable, then

Therefore, if we define $\hat{u}_n = \sqrt{n}(\hat{\beta}_{\tau_k}-\beta_0)$, then $\hat{\beta}_{\tau_k} = \beta_0+\frac{\hat{u}_n}{\sqrt{n}}$. By some simple calculations and (A.2), we have

According to condition C4, $\Sigma_1$ is a symmetric positive definite matrix. Lemma 3 is proved by Lemma 1 and Slutsky's theorem.

Proof of Theorem 4.1. By Lemma 3 we know that

From $\hat{\beta}_w = \sum_{k = 1}^{K}\omega_k\hat{\beta}_{\tau_{k} }$, $\sum_{k = 1}^{K}\omega_k = 1$, we can get

According to the proof of Lemma 2, we can obtain that

Let $\mu_i = \frac{\delta_i}{\pi(X_i, \gamma)}X_i\sum_{k = 1}^{K}\frac{\omega_k\Psi_{\tau_k}(\varepsilon_i-b_{0k})}{g(\tau_k)}, \hat{\mu}_i = \frac{\delta_i}{\pi(X_i, \hat{\gamma})}X_i\sum_{k = 1}^{K}\frac{\omega_k\Psi_{\tau_k}(\varepsilon_i-b_{0k})}{g(\tau_k)}$, and then

where $H_n = \left[\frac{1}{n}\sum_{i = 1}^{n}\frac{\partial \mu_i}{\partial \gamma} \right]_{\gamma^\ast}$. Therefore, Eq (A.4) is equivalent to

Therefore,

where $\Lambda = E[\lambda\lambda^T]$, $\lambda = \mu-E[\partial \mu/\partial \gamma^T]\{E[\partial U(\gamma)/\partial \gamma^T]\}^{-1}U(\gamma)$, $\mu = \frac{\delta}{\pi(X, \gamma)}X\sum_{k = 1}^{K}\frac{\omega_k\Psi_{\tau_k}(\varepsilon-b_{0k})}{g(\tau_k)}$.