Statistical inference for nonignorable missing-data problems: a selective review

ABSTRACT

Nonignorable missing data are frequently encountered in various settings, such as economics, sociology and biomedicine. We review statistical inference for nonignorable missing-data problems, including estimation, influence analysis and model selection. For estimation of mean functionals, we review semiparametric method and empirical likelihood (EL) approach. For estimation of parameters in exponential family nonlinear structural equation models, we introduce expectation-maximisation algorithm, Bayesian approach, and Bayesian EL method. For influence analysis, we investigate the case-deletion method and local influence analysis method from the frequentist and Bayesian viewpoints. For model selection, we present the modified Akaike information criterion and penalised method.

Keywords:

• Bayesian method
• empirical likelihood method
• expectation-maximisation algorithm
• local influence analysis
• missing data
• model selection

1. Introduction

In many applications, such as economics, sociology and biomedicine, some subjects may missing responses or predictors due to various reasons such as study drop-out, study subjects' unwilling to answer some items on a questionnaire, loss of information caused by uncontrollable factors, among others. Statistical inference for missing-data problems is quite challenging.

Over the past years, various statistical methods have been developed to handle missing data based on three different missingness data mechanisms (Little & Rubin, 2002): missing completely at random (MCAR), where missingness does not depend on any observed or unobserved data; missing at random (MAR), where missingness only depends on the observed data; missing not at random (MNAR), where missingness is associated with the unobserved data, perhaps in addition to the observed data, which is the most difficult to handle. Generally, the MCAR and MAR are called ignorable missingness, while the MNAR is termed nonignorable missingness.

The most commonly used method in missing data analysis is the complete-case analysis (i.e., statistical inference is performed only on the completely observed data). However, this method may be biased or inefficient when missingness data mechanism is not MCAR. Another popular approach to handle missing data is the imputation, including single imputation such as the hot-deck imputation (i.e., a missing value is imputed from a randomly selected similar record), the cold-deck imputation, the mean imputation (i.e., a missing value is imputed with the mean of that variable for all other observed cases), which has some attractive properties for univariate analysis but becomes problematic for multivariate analysis in that it attenuates correlation between the variable to be imputed and regression imputation (i.e., a fitted value from the considered regression model is used to impute the missing value), multiple imputation (MI) and expectation-maximisation (EM) method (Dempster, Laird, & Rubin, 1977). In particular, Horvitz & Thompson's (1952) inverse probability weighting (IPW) method is also a popular method, which may be unstable due to the usage of inverse probabilities when the probabilities are quite small even if it is unbiased under MAR, and may also suffer from efficiency loss due to incomplete data omitted. To incorporate additional information from incomplete subjects, Robins, Rotnitzky, and Zhao (1994) proposed the augmented inverse probability weighted (AIPW) method, which is also referred to as the propensity score (PS) weighting approach. Most of existing work associated with the PS method mainly focuses on the MAR assumption. When the missing data are MNAR, standard PS methods developed for ignorable missing data may not work in that missing variables in the PS model are not always observed. The PS analysis with nonignorable missing data is quite challenging. To this end, many modified PS methods have been proposed to handle missing data over the past years, for example, see Hansen (1982), Chang and Kott (2008), Kim and Shao (2013), Riddles (2013), Wang, Shao, and Kim (2014), Riddles, Kim, and Im (2016), and Jiang, Zhao, and Tang (2016).

Model-based approaches provide additional tools for handling missing data including ‘testing’ missing data types and estimating parameters in the considered models. The widely used methods include the maximum likelihood (ML) or quasi-likelihood method directly maximising the likelihood or quasi-likelihood of the observed data, and the EM algorithm for evaluating the ML estimates (MLEs) from the complete-data likelihood when the likelihood function or quasi-likelihood function is available, which has extended to various statistical models, for instance, see Laird and Ware (1982), Ibrahim (1990), Ibrahim, Chen, and Lipsitz (2001), Lee and Zhu (2002), Zhao and Tang (2015), and Tang and Tang (2018).

Also, the empirical likelihood (EL) method, which is proposed by Owen (1988), is widely applied to missing data analysis due to its attractive properties. For example, it has sampling properties similar to the bootstrap, but instead of resampling. There are considerable literature on missing data analysis, including estimations of mean functionals and model parameters using the EL method. For instance, see Cheng (1994), Wang and Rao (2002), Liang, Wang, and Carroll (2007), Stute, Xue, and Zhu (2007), Xue (2009), Qin, Zhang, and Leung (2009), Wang and Chen (2009), Tang and Zhao (2013), Wang and Qin (2010), Chen and Kim (2017), Zhao, Tang, and Tang (2013), Zhao, Zhao, and Tang (2013), and Tang, Zhao, and Zhu (2014).

Bayesian parameter estimation method is also often used to handle the estimation problem of the parametric models with missing data due to the following facts. First, the posterior distributions of parameters and nonparametric function and missing data can be estimated using the well-known Markov chain Monte Carlo (MCMC) methods in statistical computations, such as the Gibbs sampler (Geman & Geman, 1984) and the Metropolis-Hastings (MH) algorithm (Hastings, 1970; Metropolis, Rosenbluth, Rosenbluth, Teller, & Teller, 1953). Second, Bayesian method with missing data only requires an additional step in the Gibbs sampler compared with no missing-data settings. Thus, Bayesian method can easily accommodate missing data without requiring new techniques for statistical inference (Ibrahim, Chen, Lipsitz, & Herring, 2005). Third, some prior information can be directly incorporated into the analysis, which yields more precise estimates of parameters in settings with good prior assumption. Fourth, the sampling-based Bayesian methods do not depend on asymptotic theory, and may give more reliable statistical inference even with small sample sizes (Lee & Tang, 2006a). There are a lot of discussions on Bayesian analysis for missing data problems. For example, see Ibrahim et al. (2005), Lee and Tang (2006b), Lee and Tang (2006c), Daniels, Wang, and Marcus (2014), Linero and Daniels (2015), Tang, Chow, Ibrahim, and Zhu (2017), Zhang and Tang (2017).

The aforementioned model-based methods heavily depend on the missingness data mechanism and modelling assumptions, thus the resulting estimators may be sensitive to these assumptions. Hence, sensitivity analysis for the missing data problems, which is usually used to investigate the effect of the data and model assumptions on estimation outputs, has received considerable attention over the past years. To address the issue, there are two major approaches: case-deletion method (Cook, 1977) and local influence analysis (Cook, 1986), which have been extended to various statistical models in the absence of missing data (e.g., see Cook & Weisberg, 1982; Wei, 1998). Also, the two approaches have been extended to the missing data problems. For example, see Lavine (1991), Berger (1994), Dey, Ghosh, and Lou (1996), Gustafson (1996a,1996b), Zhu and Lee (2001), Verbeke, Molenberghs, Thijs, Lesaffre, and Kenward (2001), Jansen, Molenberghs, Aerts, Thijs, and Van Steen (2003), Troxel, Ma, and Heitjan (2004), Millar and Stewart (2007), Daniels and Hogan (2008), Shi, Zhu, and Ibrahim (2009), Zhu, Ibrahim, Cho, and Tang (2012), Zhu, Ibrahim, and Chen (2015). In particular, Zhu, Ibrahim, and Tang (2011) developed a unified approach of Bayesian influence analysis for simultaneously assessing the effect of various perturbations to the data, the prior and the sampling distribution for a class of statistical models. Several extensions of Zhu, Ibrahim, and Tang (2011) have been considered in recent years, for example, see Kaciroti and Raghunathan (2014), Zhu, Ibrahim, and Tang (2014), and Zhang and Tang (2017).

However, in practical applications, the assumed missing data mechanism models are not testable from the data due to the missing value involved. To this end, various approaches have been developed to select an appropriate model, which fits the data best among the candidate models. This is a growing, active and open research field. The commonly used approaches to select the best fitting model include the stepwise regression and best subset selection associated with the Akaike's information criterion (AIC) (Akaike, 1973) and Bayesian information criterion (BIC) (Schwarz, 1978) in the frequentist framework, which heavily depend on the observed data likelihood. In the presence of missing data, it is rather difficult to obtain the observed data likelihood due to the intractable multiple integration involved. In this case, many model selection methods have been developed for missing data based on the EM algorithm or imputation technique. For example, see Hens, Aerts, and Molenberghs (2006), Ibrahim, Zhu, and Tang (2008), Shen and Chen (2012, 2013), Long and Johnson (2015), and Jiang, Nguyen, and Rao (2015).

When the number of predictors is large, the aforementioned model selection strategies often suffer from instability and heavily computational problem in the presence of missing data. To overcome the difficulties, various regularisation methods have been proposed to simultaneously estimate parameters and select important predictors over the past years. For example, see Fang and Shao (2016), and Tang and Tang (2018).

The rest of this paper is organised as follows. Section 2 introduces missing data mechanisms. Section 3 reviews two approaches to estimate the mean functionals in the presence of nonignorably missing response data. Section 4 investigates the parameter estimation problem in exponential family nonlinear structural equation models (EFNSEMs) with nonignorable missing data including EM algorithm, Bayesian method and Bayesian EL approach. Section 5 introduces influence analysis methods in the presence of missing response not at random. Section 6 discusses the variable selection problem in the presence of missing response not at random. Some concluding remarks are given in Section 7.

2. Missing data mechanism

In incomplete data analysis, the first problem being faced is how to deal with missing data. Most works in early missing studies suggest using the complete and available cases, i.e., omitting/deleting these incompletely observed subjects, to make statistical inference. This opinion essentially thinks that subjects with missing values are just a random subset of the complete sample of subjects. However, it is unreasonable when the data present the dependent structure. This issue has not been attracted researchers' attention until Rubin's pioneering work Rubin (1976). Rubin (1976) suggested that each incomplete data has its own probability generating mechanism and each observation according to some particular probability given the incomplete data is missing or observed. This opinion naturally admits that the ‘true’ values of complete data exist, but some values are observed while the others are missing with some uncertain reasons. Thus, modelling such uncertainty is a key work in missing data analysis. To efficiently quantify such uncertainty of missing value, Rubin (1976) first defined three missing data mechanisms: missing completely at random (MCAR), missing at random (MAR) and missing not at random (MNAR). In accordance with Little and Rubin (2002), missing data are ‘MCAR’ if the failure to observe a value does not depend on any observed or missing values; missing data are ‘MAR’ if the failure to observe a value only depends on the observed values rather than any unobserved values; missing data are ‘MNAR’ if the failure to observe a value depends on the value that would have been observed. This modelling strategy does not intend to explore the reasons leading to missingness (in fact, accurately describing all the potential causes or reasons for missingness is not realistic). Instead, it adopts a mathematical device to describe the rates and patterns of missing values and to roughly capture possible relationships between the missingness and the values of missing variables themselves. In this sense, it is consistent to statistical sayings ‘learning from what the data speak’. To comprehensively understand the definition of missing data mechanism, a formal mathematical form of missing data mechanism is given as follows.

Let X be a vector of explanatory variables and Y be a response variable, and be n independent realisations of continuous random variables from a joint distribution . Here, we assume that is always observed but is subject to missingness. Let be the missing data indicator for , that is, if is missing, and otherwise. Throughout this paper, we assume that is a Bernoulli random variable with the respondent probability , where is referred to as the propensity score of missing response, and is independent of for any . Thus, we assume that the missing data mechanism for is given by When 's are unrelated with the data , the missing data mechanism is called ‘MCAR’, which is the simplest missing data mechanism. For the MCAR mechanism, subjects who have missing data are just a random subset of the complete sample of subjects, and the widely used technique for handling missing data is the complete-case analysis yielding unbiased conclusions. When 's only depend on rather than (i.e., ), the missing data mechanism is termed ‘MAR’, which is a commonly used assumption in missing data literature. Under the MAR assumption, the complete-case analysis yields biased results. Thus, some new methods have been developed to handle missing data, for example, see Cheng (1994) and Wang and Rao (2002). When 's depend on perhaps in addition to , the missing data mechanism is referred to as MNAR, which is also called nonignorable. MNAR data analysis is rather challenging due to missing value involved. Using MAR data analysis method to make inference on MNAR data may lead to biased results. To address the issue, some new methods for handling MNAR data have been proposed in recent years, for example, see Qin, Leung, and Shao (2002), Kim and Yu (2011) and Tang et al. (2014).

Under some mild conditions, Rubin (1976) showed that statistical inferences based on the observed data and missing information specified by MCAR and/or MAR are the same as those based on the observed data, and hence the MCAR and MAR missing data mechanisms are called ignorable. Under MNAR missing data mechanism, efficient inferences usually require correctly specifying either missing data mechanism model or distribution of , or both. Thus, resulting conclusions may be sensitive to these assumptions. Therefore, correctly specifying missing data mechanism is rather important in MNAR data analysis but is quite difficult in practical applications because of the missing data mechanism is not testable.

3. Estimation of mean functionals with MNAR

In this section, our main interest is to estimate in the presence of missing response not at random.

3.1. Semiparametric estimation

In the nonignorable missing data literature, the missing data indicator is usually formulated by a fully parametric model, such as logistic model, probit model and log-linear model, based on the assumption that covariates has a linear effect on the function of the respondent probability for unobserved . In practice, this assumption may be inappropriate in some settings because the true relationship between the respondent probability for and covariate is usually unknown. Moreover, the fully parametric model is sensitive to the misspecification of the posited model, and is not testable under nonignorable missingness (Ibrahim et al., 2001). To solve the aforementioned problem, Kim and Yu (2011) defined the following exponential tilting model for the propensity score: (1) where is some unknown function, and φ is an unknown titling parameter. Thus, (1(1) ) defines a semiparametric model and is weaker than the fully parametric model given in Lee and Tang (2006b). To estimate the propensity , Kim and Yu (2011) considered two cases that φ is known or can be estimated using the additional data. When , model (1(1) ) reduces to an MAR model. Let be the conditional density of given and , and let be the conditional density of given and . Then, following the argument of Kim and Yu (2011), we have (2) where the tilting parameter measures the amount of departure from the MAR assumption. When , (2(2) ) reduces to , which shows that our considered missingness data mechanism is a natural extension of MAR mechanism, and MAR mechanism is a special case of our considered missingness data mechanism.

Under the assumption of model (1(1) ), and γ are unidentifiable when both are unknown Shao and Wang (2016). To address the identifiability with an unknown tilting parameter φ, Shao and Wang (2016) show that all unknown parameters can be identified and consistently estimated if one can find an instrument (e.g., z) that is not emerged in (1(1) ), i.e., if and (3) where is still an unknown function of u.

Following Shao and Wang (2016), when γ is known, we can rewrite (1(1) ) as . Thus, we have (4) which leads to When γ is the true parameter, it follows from Kim and Yu (2011) that the kernel regression estimator for can be expressed as where , and is a symmetric kernel function on the real line and is a smoothing bandwidth sequence such that and as .

When γ is unknown, under model (3(3) ), Equation (4(4) ) still holds with x replaced by u. It is assumed that the instrument variable z takes the finite discrete values such as . For each l, we define where is the indicator function. Denote , . Then, we have (5) which is called instrumental estimating equations Shao and Wang (2016). Note that when L=1, there is no instrumental variable, Equation (5(5) ) reduces to (4(4) ) and is not sufficient for estimating and γ.

To estimate , the idea of profiling is here adopted. For every fixed γ, the sum of L equations in (5(5) ) is used to obtain the kernel estimator of , which is given by where the bandwidth depends on l. Note that is not an estimator of unless γ is the true parameter, but it is useful for profiling. More importantly, estimating only uses the sum of L equations, the left L−1 estimating equations are used to estimate unknown parameter γ Shao and Wang (2016), which shows that we require . After we obtain , substituting for in (5(5) ) leads to the following profiled L estimating equations: (6) A sample version of the profiled instrumental estimating Equation (6(6) ) is (7) When L is greater than , where is the dimension of γ, (7(7) ) has no solution due to its over-identifying. In this case, a two-step generalised method of moments Hansen (1982) is adopted to estimate γ. To save space, we omit it. The details see Shao and Wang (2016).

When the missing data mechanism for Y is not MAR but MNAR, the estimator with Cheng (1994) is biased. In this case, θ can be estimated by , where is a consistent estimator of . Under MNAR assumption, estimating is quite challenging in that is not obtained in the set of nonrespondents. When γ is known, Kim and Yu (2011) proposed a nonparametric regression estimator for , which is given by where the weight presents the point mass assigned to when is approximated by . Thus, a nonparametric estimator of θ can be obtained by When , estimator reduces to . Under some regularity conditions, Kim and Yu (2011) showed where . The increase in variance due to missing data is , where . The above equation indicates that the variance increase depends on the respondent probability and the squared error . When , the variance increase is zero. When is quite small, the variance increase can be quite large. When does not depend on Y (i.e., ), reduces to the above defined . In practical applications, is usually unknown. To this end, Kim and Yu (2011) presented a consistent estimator of , which is given by where , and with .

When γ is unknown, γ can be estimated from an independent survey or a validation sample, which is a subsample of the nonrespondents. Let be the corresponding estimator of γ. Thus, a semiparametric estimator of is given by , which indicates that a semiparametric estimator of θ can be expressed by . When is obtained from an independent survey, under some assumptions such that and is independent of , Kim and Yu (2011) showed where , and . The above equation shows that the increase in variance is due to estimating γ. A consistent estimator of can be expressed as , where is a consistent estimator of , and with . When is obtained from a validation sample, under some regularity conditions, Kim and Yu (2011) proved where , in which r is an indicator function taking 1 if subject belongs to the follow-up sample and 0 otherwise, and , and is the probability limit of . It is easily seen from the definition of that when , and is minimised due to the fact that when the assumed respondent model (1(1) ) is true. Thus, the validity of the proposed estimator does not depend on the assumed respondent model and the role of the respondent model is to improve the efficiency. A consistent estimator of is given by where .

3.2. EL estimation

Due to some good properties of EL estimation, Qin et al. (2002) developed a semiparametric EL approach to estimate mean functions for the data with nonignorable nonresponse. It is rather difficult to evaluate semiparametric EL estimate of θ due to simultaneously solving multiple equations. To address the issue, Zhao et al. (2013) developed two EL methods to estimate under MNAR assumption. In what follows, we only consider the case that γ is known. When γ is unknown, we can substitute a consistent estimator of γ to get a semiparametric estimation. Let . The estimated empirical log-likelihood ratio function for parameter θ is given by When with probability tending to one, the optimal value of is , where is the Lagrange multiplier and satisfies . Thus, we have The maximum EL estimator (denoted by ) of θ can be obtained by maximising , which shows that under some smooth conditions, can be obtained by simultaneously solving equations: and . The nested optimisation routines (Owen, 2001, Ch. 12) can be here adopted to evaluate . When , reduces to the EL estimator of Wang and Rao (2002).

Under MNAR assumption, it is easily shown that . Hence, the aforementioned maximum EL estimator for θ based on may be biased. To address the issue, Zhao et al. (2013) proposed the augmented IPW imputation procedure. Let . Under MNAR assumption (1(1) ), we have , which shows that the random variable is an unbiased estimator of θ. In applications, is usually unknown because of a nonparametric function involved. Similar to Kim and Yu (2011), we can substitute for to get a semiparametric estimator of as . Then, the estimated empirical log-likelihood ratio function for θ is defined as By using the Lagrange multiplier method, can be expressed as , where satisfies . Thus, the maximum EL estimator (denoted by ) of θ can be obtained by maximising .

It is assumed that some auxiliary information on X may be available. For instance, the mean of X is equal to zero or the distribution of X is symmetric. In particular, we assume that the auxiliary information on X has the form of , where is a known vector function. In this case, the estimated empirical log-likelihood ratio function for θ can be defined as where . Using the Lagrange multiplier method, we obtain where is the solution to the following equation: . Similarly, the maximum EL estimator (denoted by ) of θ based on the auxiliary information can be obtained by maximising.

Under some regularity conditions, if θ is the true parameter, Zhao et al. (2013) showed

where , , , and with and .

The above equations characterise the asymptotic normality of , and , which can be employed to construct the normal-approximation-based (NA-based) confidence interval or region for θ. For instance, an approximate 100( NA-based confidence region for θ is given by . They also indicate that and have the same asymptotic variance, which can be explained by the fact that the nonparametric kernel estimator of the response probability is a consistent estimator when γ is known. Based on the asymptotic normality of maximum EL estimators, we can approximate covariance matrices of , and . Similarly, a consistent estimator of is given by . The consistent estimators of and are given by and , respectively. Thus, the consistent estimators of and are , and , respectively.

Also, it is easily known from the above equations that (i) has smaller asymptotic variance and hence is asymptotically more efficient than and ; (ii) when the auxiliary information on X is available and X is not subject to missingness, the reduction in asymptotic variance of and and that of does not depend on the respondent probability function , which is consistent with that under MAR assumption for Y (Wang & Rao, 2002). When the missingness data mechanism for Y is MAR, we have , which indicates that and reduce and , respectively. Thus, Theorems 2.1 and 3.1 in Wang and Rao (2002) and Theorems 1–3 in Xue (2009) have been extended to our considered case.

We can use the asymptotic distributions of , and to construct the confidence region of θ. For example, let be the quantile of for . Thus, an approximate 100 EL-based confidence region for θ is defined by CI.

4. Estimation of parameters in EFNSEMs with MNAR

In educational, behavioural, social and psychological studies, structural equation models (SEMs) are widely used to identify relationships among observed and latent variables (Bentler & Wu, 2002; Jöreskog & Sörbom, 1996; Lee & Tang, 2004; Lee & Zhu, 2002). Many structural equation modelling packages have been developed to handle the strong demand in various fields; for example, see LISREL (Jöreskog & Sörbom, 1996), EQS6 (Bentler & Wu, 2002), Mplus (Muthén & Muthén, 2017), WinBUGS (Spieglhalter, Thomas, Best, & Lunn, 2003). These packages mainly focus on normal data. In many applications, there are a large number of non-normal manifest variables. To this end, it is important to develop statistical inference tools to analyze NSEMs with manifest variables following exponential family distribution, which includes many commonly encountered distributions such as the normal, exponential, binomial and Poisson distributions as its special cases.

4.1. MLE

For , let be a vector of manifest variables measured on each of n subjects, and be a vector of latent variables. To measure the relationship between and , it is assumed that components of given are conditionally independent, and () given and is distributed as the following one-parameter exponential family (8) with and , where and are specific differentiable functions. Exponential family distribution includes many distributions as its special cases, such as the normal, exponential, binomial and Poisson. Similar to generalised linear models (McCullagh & Nelder, 1989), we measure the relationship between and via the following measurement model: (9) where u is a vector of unknown parameters, is a vector of explanatory variables, β is a matrix of regression coefficients, Λ is a factor loading matrix. To identify the relationships among the latent variables, we partition as and assume that and satisfy the following nonlinear structural equation: (10) where and are and latent subvectors of , respectively, and ; ϒ is a matrix of regression coefficients, is a vector of covariates, is a nonzero vector-valued function with differentiable functions and ; Π () and Γ () are matrices of unknown regression coefficients of on and . Further, it is assumed that is a nonzero constant and independent of Π, and are independently distributed as and , respectively; where . Let and . Thus, the nonlinear structural Equation (10(10) ) can be rewritten as.

The above defined models (or their parameters) are not identifiable, and some identification conditions on these parameters are required. Following the widely adopted technique in SEM (Lee & Tang, 2006b), we fix some components of Λ and/or to some known values, for example, 0 or 1.

For simplicity, we assume that the manifest vector is incompletely observed and subject to nonignorable missingness, and and are fully observed. We define a missing indicator , which takes 0 if is missing, and 1 if is observed. Similar to Lee and Tang (2006b), it is assumed that the conditional distribution of and given and are independent for . Under this assumption, we consider the following model for nonignorable missingness mechanism: (11) where , , and and ϕ is a appropriate parameter vector. Since is binary, a sequence of logistic regressions can be adopted to model . Similar to Lee and Tang (2006b), we consider the following logistic regression model for : (12) (LoG) where , and and . The missingness data mechanism defined in (12(12) (LoG) ) is called missing not at random(MNAR). More discussions on the specification of missingness data mechanism models can see Lee and Tang (2006b).

Denote , , , , , , , and . Let θ be a parameter vector containing all unknown distinct parameters in , , Φ, and ϕ. The joint probability density function of θ based on the complete data is given by where , is the probability density function of given , and is the probability density function of given . Denote . Thus, the log-likelihood function of θ based on the observed data is given by which is rather complicated and has not closed-form. Hence, it is quite difficult to compute the ML estimation of θ based on . To address this issue, the well-known EM algorithm is adopted to evaluate the ML estimation of θ.

Based on the idea of data augmentation, the observed data is augmented with the missing data and latent factors ω. Thus, we obtain the complete data set . Based on the missingness data mechanism model (12(12) (LoG) ), the complete-data log-likelihood function for θ is given by (13) (LogLK) where (14) with , (15) (LogLK2) (16) (LogLK3) (17) (LogLK4) where , and represent the constant that don't depend on parameters. Clearly, is simpler than in that it does not involve intractable multiple integrals, and has four separate terms with separate parameters. Let be the conditional expectation of with respect to the conditional distribution , where represents the MLE of θ. The EM algorithm for evaluating the MLE of θ usually requires repeating the following two steps: E-step (i.e., Expectation step) and M-Step (i.e., Maximum step) until the algorithm converges. That is, for at the κth iteration of the EM algorithm with a current value of θ (e.g., corresponds to the initial value of θ),

E-step: compute , where the expectation is taken with respect to the conditional distribution of given ;

M-step: find by maximising with respect to θ.

It is easily seen from the definition of that it is still quite difficult to compute because of the multiple integrals involved. To this end, the well-known Monte Carlo method is employed to approximate . Following Wei and Tanner (1990), Monte Carlo approximation to can be expressed as with , where is a sample drawn from the conditional distribution  via the Metropolis-Hastings (MH) algorithm within the Gibbs sampler. Here, the Gibbs sampler is implemented as follows. Given the value of θ, at the tth iteration of Gibbs sampler with the current values and of and ω, we first generate from the conditional distribution , and then draw from the conditional distribution .

In what follows, we derive the conditional distributions and . Since 's are conditionally independent, and 's are also conditionally independent and do not depend on Z, , Φ and , it follows from (8(8) ) and (9(9) ) that (18) (CondYmis) which indicates that the conditional distribution is proportional to (19) with , where . The MH algorithm for drawing observations from the conditional distribution is as follows. At the tth iteration of Gibbs sampler with the current values and of and , a candidate is generated from the normal distribution , where is selected with the average acceptance rate being roughly 0.25, and with being the subvector of corresponding to missing components of , , being the subvector of ϕ corresponding to , being the subvector of ϕ corresponding to , and .

It follows from the independence of subjects that . It follows from the assumptions given in Equations (8(8) ), (9(9) ) and (12(12) (LoG) ) that (20) (Condomga) with . Similarly, the MH algorithm for sampling observations from is implemented as follows. At the tth iteration of Gibbs sampler with the current values and of and , a candidate is generated from the normal distribution , where is selected as , and with and .

In the M-step, we need maximising with respect to θ. However, there is no analytic solution of θ to . To this end, the widely used ECME algorithm (Liu & Rubin, 1994) is used to separately maximise four terms in at every conditional maximum (CM) step. Let be the kth row vector of matrix , and be the kth component of vector . Thus, the closed forms for , and Φ at the th conditional maximum likelihood (CML) step of the EM algorithm with the current values , and are given by respectively.

But, there are no analytic solutions for u, , and ϕ at the κth CML step. To this end, Newton-Raphson maximisation procedure is adopted to evaluate the solutions of the intractable CM step for u, , ψ and ϕ because of its attractive selection in the local optimisation. Let be the jth row subvector of . Define , , and . Denote , , , , where Given the current estimates , , and of u, , and ϕ at the κth iteration of the EM algorithm, Newton-Raphson iterative formulae for evaluating , , and are given by respectively, where , , , and .

When Newton-Raphson algorithm converges, the corresponding value of is denoted as , which is the desirable ML estimates of θ. The above Newton-Raphson iterative procedure is repeated until the infinite norm is sufficiently small, for example, 0.001. Convergence of the above presented ECME algorithm is monitored by using the log-likelihood ratio statistic via the bridge sampling method (Meng & Schilling, 1996).

Similar to missing data literature (e.g., see Lee & Zhu, 2002), we can use the inverse of the observed information matrix to approximate the asymptotic covariance matrix of . However, it is rather difficult to obtain the analytic expression of the observed information matrix because of the multiple integral involved. Here, an identity of Louis (1982) is adopted to evaluate standard error estimates of components in on the basis of the generated observations from the conditional distribution via the above presented MH algorithm within the Gibbs sampler. Following Louis (1982), can be approximated by where , and for any vector a. Expressions of and are given in the Appendix.

4.2. Bayesian estimation

According to the above argument, it is rather difficult to make Bayesian inference on θ via the observed data likelihood function . Inspired by the idea of data augmentation (Tanner & Wong, 1987), we regard as hypothetical missing data, and augment them with the observed data set in the posterior analysis. The posterior distribution , which is proportional to , is simpler to evaluate than , which is the posterior distribution based on the observed data, where is the prior density of θ. However, it is rather difficult to make inference on θ from the posterior distribution because it is not easy to directly draw observations from . To address the issue, Bayesian analysis is conducted on the basis of the complete data density function through the MCMC methods, such as the Gibbs sampler (Geman & Geman, 1984) and the MH algorihtm. Here, the Gibbs sampler is adopted to draw a sequence of random observations from the posterior distribution , which is conducted as follows. At the th iteration of the Gibbs sampler with current values , , and , we iteratively

1. draw observation from the conditional distribution ,

2. draw observation from the conditional distribution ,

3. draw observation from the conditional distribution . Following the argument of Geman and Geman (1984), under some regularity conditions and for a sufficiently large , observations can be regarded as an observation from the posterior distribution . Hence, the posterior distribution can be approximated by the empirical distribution of observations , where is a sufficiently large number. The widely used ‘estimated potential scale reduction’ (EPSR, Gelman, 1996) values of the parameters in θ, which are computed sequentially as the iterations continue from several different starting values of parameters, is employed to monitor the convergence of the Gibbs sampler. When the EPSR values of parameters in θ are less than 1.2, we claim that the algorithm converges.

(1) Prior specifications

For the parameters u, β and Λ given in the measurement model (9(9) ), we consider the following priors: (21) where and () are hyperparameters whose values are prespecified, and Gamma represents the Gamma distribution with parameters a and b. For , it is assumed that is independent of . When and are selected to be large, and is taken to be small (e.g., 0.01), the above given priors reduce to the situation with little information.

For the parameters , and Φ associated with the structure equation model (10(10) ), we consider the following prior: (22) where , , , , and are hyperparameters whose value is assumed to be pre-given, and denotes the inverted Wishart distribution. When approaches to 0, it reduces to a non-informative prior.

For the parameter vector ϕ associated with the missing mechanism model (12(12) (LoG) ), a normal prior is specified: (23) where and are hyperparameters. Similarly, when approaches to 0, it reduces to a non-informative prior.

(2) Conditional distributions

Under the aforementioned assumptions, the joint distribution of has the following decomposition: It is reasonable to assume from the definition of θ that where , , , , , and are the prior distributions of and ϕ, respectively. Thus, the joint distribution of can be written as (24) (Joint) which is adopted to derive the conditional distributions requiring in implementing the Gibbs sampler as follows.

For the conditional distribution given in (18(18) (CondYmis) ), it follows from independence of subjects that . While the conditional distribution is given in (19(19) ).

For the conditional distribution given in (20(20) (Condomga) ), since is independent of for , it follows from (24(24) (Joint) ) that . While the conditional distribution is given in (20(20) (Condomga) ).

Consider the conditional distributions , and . It follows from (8(8) ), (9(9) ), (24(24) (Joint) ) and (21(21) ) that (25) (Postu) (26) (PostBeT) (27) (PostPSI) for , where and . Particularly, if does not depend on , where is some appropriate function, we have .

Consider the conditional distributions , and . It follows from (10(10) ), (22(22) ) and (24(24) (Joint) ) that the joint conditional distribution is proportional to which leads to (28) where , and . Also, it is easily shown from (10(10) ), (22(22) ) and (24(24) (Joint) ) that the conditional distribution is proportional to which yields (29) Consider the conditional distribution . It follows from (12(12) (LoG) ), (23(23) ) and (24(24) (Joint) ) that (30) Note that the conditional distributions of and are given under the situation without fixed parameters. The method presented in Lee and Tang (2006b) can be applied to handle the situation with fixed parameters.

(3) Implementation

It is easily seen from (28(28) ) and (29(29) ) that conditional distributions , and are the standard distributions. Thus, it is fast and easy to sample observations from these conditional distributions. But, conditional distributions , , , , , and given in (19(19) ), (20(20) (Condomga) ), (25(25) (Postu) )–(27(27) (PostPSI) ) and (30(30) ), respectively, are non-standard and complicated distributions. Hence, it is rather difficult to directly draw observations from these conditional distribution. To address the issue, the MH algorithm is adopted to generate observations from these conditional distributions. The MH algorithms for sampling and from their corresponding conditional distributions have been given in Section 4.1.

To sample observation u from via the MH algorithm, we take as the target density, and select as the proposal distribution, where is selected so that the average acceptance rate is approximately 0.25, and with in which and for . The MH algorithm is implemented as follows. At the tth iteration of the Gibbs sampler with a current value , a new candidate u is sampled from and accepted with probability Similarly, the MH algorithms for sampling observations from , and are as follows. At the tth iteration of the Gibbs sampler with the current values , and , new candidate values , and ϕ are generated from the proposal distributions , and , respectively, and accepted with probabilities respectively, where and , , and are chosen so that the average acceptance rate is approximately 0.25.

(4) Bayesian estimates of parameters and latent variables

Based on observations sampled from the aforementioned presented MH algorithm together with the Gibbs sampler procedure, we can obtain the joint Bayesian estimates of parameters in θ and latent variables 's and their corresponding standard errors.

Let be the sampled observations of from their joint conditional distribution via the aforementioned presented MCMC algorithm. Thus, Bayesian estimates of and can be evaluated by respectively. Similarly, the consistent estimate of the posterior covariance matrix var can be computed by the sample covariance matrix of observations drawn using the aforementioned presented MCMC algorithm. For example, . Hence, the estimated standard errors for components of θ can be evaluated from the diagonal elements of the sample covariance matrix.

4.3. Bayesian EL estimation

Although the aforementioned given assumptions are common and reasonable in many cases, it may lead to biased estimates or even misleading conclusions when the true distribution of response variable is highly non-exponential. In these cases, we assume that the distribution of is unknown, but satisfying , and . It is also assumed that 's are subject to non-ignorable missingness, and 's are completely observed. The assumption for is the same as that given in Section 4.1.

Let . When 's are completely observed, under the above assumption, we can construct the following estimating equations: which satisfies , where . Under the assumption given in Equation (10(10) ), we can construct the following estimating equations: which satisfies . Under the assumption , we can construct the following estimating equations: which satisfies . According to the assumption given in Equation (12(12) (LoG) ), we can construct the following estimating equations: which satisfies , where . Also, under the assumption given in Equation (10(10) ), we can construct the following estimating equations: which satisfies .

Denote , , , . Following Owen (2001), we define the profile EL ratio function for θ as It follows from Owen (2001) that the optimal value of can be expressed as (31) where is the root of the following equation: Thus, we have . Hence, the EL function of θ is given by , which can be regarded as the working likelihood of θ. Based on the working likelihood of θ and the prior distribution of θ, the posterior probability density of θ can be expressed as , where .

Under our considered situation that 's are subject to missingness and 's are unobserved latent variables, it is impossible to make Bayesian inference on θ based on . To address the issue, a hybrid algorithm is proposed to impute latent variables and missing responses by combining the Gibbs sampler and the MH algorithm based on the posterior distributions , and . In the Gibbs sampler, the first difficulty is the selection of the initial values of unobserved latent variables and missing responses. A simple approach for solving the issue is to independently generate several MCMC observations based on the classical SEM with the complete-case analysis. The second difficulty is to estimate the conditional distributions and . To this end, the approach proposed by Wei, Ma, and Carroll (2012) is employed to estimate and based on the EL , and then latent variables and missing responses are drawn by using the linear interpolation method. Following Zhang and Tang (2017), the algorithm for implementing the Gibbs sampler can be summarised as follows.

Step 1. Let t=0 and take the initial values and of θ and , respectively.

Step 2. Sample observations of θ from its posterior distribution. The MH algorithm for sampling θ is conducted as follows. At the th iteration with a current value of θ, we first draw a candidate value of θ from the proposal distribution , where is chosen such that the average acceptance rate is roughly 0.3; and then accept as with the probability where 's are generated by using the modified Newton-Raphson algorithm (Chen, Sitter, & Wu, 2002). Note that it is impossible to implement the above procedure because components of θ have different features. In this case, we can iteratively implement the above procedure for u, , , Φ, ϕ and .

Step 3. Sample latent variables from the corresponding conditional distributions. First, we estimate the conditional density . To this end, we take grid points (e.g., ) in some neighbourhood of , and calculate for based on Equation (31(31) ). The conditional density function of latent variable is estimated by where , is the ιth component of . Thus, the cumulative distribution function of latent variable can be estimated by Following the idea of the linear interpolation method (Wei et al., 2012), we first sample a random variable u from the uniform distribution , and then take as the th iteration value of . Thus, .

Step 4. Sample missing response from the corresponding conditional distribution. First, we estimate the conditional density . To this end, we take grid points (e.g., ) in some neighbourhood of , and calculate for based on Equation (31(31) ). The conditional density function of is estimated by where , is the ιth component of . The cumulative distribution function of missing response is estimated by for . Similarly, we first sample a random variable u from , and define as the th iteration value of . Thus, .

Step 5. Repeat Steps 2–4 until convergence of the algorithm. Bayesian estimates and of θ and can be taken to be their corresponding sample means of the generated observations.

5. Influence analysis of EFNSEMs with MNAR

5.1. Local influence analysis

Let be a vector of perturbations restricted to some open subspace Ξ, and be the complete-data log-likelihood of the perturbed model. Let be the ML estimate of θ for the perturbed model. That is, is obtained by maximising . We assume that there is a such that for all θ. Thus, we have . To wit, represents no perturbation. Following Zhu and Lee (2001), we consider the following Q-displacement function . Similar to Cook (1986), we define the influence graph of as . Thus, the normal curvature of in the unit direction () at can be defined as where , , and . It follows from Poon and Poon (1999) that the conformal normal curvature in the unit direction at can be expressed as (32) Let , and be the nonzero eigenvalues of , and be the corresponding orthogonal eigenvectors. Similar to Cook (1986), the eigenvector corresponding to the largest eigenvalue is important for assessing local influence. However, as many authors (e.g., Lesaffre & Verbeke, 1998; Poon & Poon, 1999) pointed out that it is not enough to assess local influence by only examining . To address the issue, Poon and Poon (1999) proposed using statistic to assess local influence, where for . It is easily seen that for , where is the jth component of , is an vector of perturbations with its jth entry 1 and zeros elsewhere, and is the jth diagonal element of matrix . Thus, evaluating is rather simple because of without computing eigenvectors or eigenvalues of matrix . Let and SM be mean and standard error of , respectively. Similar to Lee and Tang (2004), is used as a benchmark, where c is a selected constant on a problem-by-problem basis. To wit, the jth observation is identified as influential if .

As Lee and Tang (2004) pointed out, it is quite difficult to evaluate and because of multiple integrals involved. To this end, Lee and Tang (2004) employed the well-known Monte Carlo method to approximate and . Let be a set of observations generated from the conditional distributions and given in Equations (19(19) ) and (20(20) (Condomga) ) via the MH algorithm given in Section 4.2.

Thus, and can be approximated by According to the above introduced EM algorithm, at the E-step of the last MCECM iteration, a sample has been generated at the ML estimate . If , no additional observations are required; otherwise the MH algorithm is used again to sample additional observations to make up a sufficient large sample for local influence analysis. Generally, is selected on a problem-to-problem basis. In some applications, a few thousand (e.g., 5000) observations are sufficient for models with moderate sample sizes.

The above presented local influence measures are mainly defined on the basis of and , which are closely associated with , and , . It follows from Equation (13(13) (LogLK) ) that the complete-data log-likelihood function , is a relatively simple function of three separate terms with separable parameters. Hence, is a diagonal block matrix, and many components of are zero for many perturbed models. Thus, the computation of the above defined local influence diagnostic measures is quite simple. As an illustration, six perturbation schemes are considered as follows.

(1) Case weights perturbation

We first consider simultaneous changes in the weights of all subjects. Let be an vector of case-weights. In this case, represents no perturbation, and the complete-data log-likelihood function for the perturbed model is written as (33) where , , . In particular, when and for , we obtain , where is the ML estimate of θ after deleting the ith data point. Thus, this perturbation scheme is an extension of the case-deletion method.

(2) Additive perturbation on components of the manifest vectors

Here we consider an additive perturbation on the manifest vector, i.e., , where are the kth component of . In this case, , , and corresponds to a case with no perturbation. The complete-data log-likelihood function for the perturbed model is given by (33(33) ), but , with . In particular, we consider a perturbation on the manifest vector , i.e., , where is a scalar, and .

For the above two perturbation schemes, and , in are nonzero, while other two terms are zero. Moreover, only involves , and and and . Thus, the computational burden for evaluating diagnostic measures is not too heavy. Also, the two perturbations give us more precise and direct insights with clear interpretations.

(3) Perturbation on

Consider perturbation on using a vector such that . In this case, represents no perturbation. The perturbed complete-data log-likelihood function is given by (34) where . It is easily seen from Equation (34(34) ) that is the unique nonzero derivative in , and computation of does not involve , and . Hence, local influence diagnostics for this perturbation scheme only depend on and , and the nonlinear SEM.

(4) Multiplicative perturbation on the latent vector

This perturbation scheme has been considered by Lee and Tang (2004). Following Lee and Tang (2004), we consider perturbation on : for , where . This perturbation scheme includes either the multiplicative perturbation on the latent subvector or by using or , respectively. In this case, represents no perturbation. The complete-data log-likelihood function for the perturbed model is given by where

(5) Additive perturbation on the latent vector

This perturbation scheme has been considered by Lee and Tang (2004). Following Lee and Tang (2004), we consider an additive perturbation on the latent vector by taking . In this case, represents no perturbation. The perturbed complete-data log-likelihood is given by where The above three perturbation schemes are closely associated with the nonlinear SEM (10(10) ) that measures the relationships among the latent variables. These perturbation schemes can be used to detect the effect of the latent variables on the manifest observations.

(6) Perturbation on all unknown parameters

Similar to Lee and Tang (2004), we consider perturbation on all unknown parameters. In this case, the perturbed parameter vectors are given by , , , , where , and are diagonal matrices of appropriate orders, and is a matrix that can be determined by users; and is the perturbed vector corresponding to for i=1,2, and and are the perturbed vectors corresponding to Φ and ϕ, respectively. In this case, represent no perturbation for i=1,2,3,4. The complete-data log-likelihood function of the perturbed model has the form of Therefore, we have where is a diagonal block matrix with blocks , , and .

5.2. Case-deletion analysis

In what follows, a quantity with a subscript ‘’ denotes the quantity evaluated from a reduced data set, which is obtained by deleting the ith observation from the full sample . For example, denotes the ML estimate of θ with the ith data point deleted. That is, is obtained by maximising the objective function , where represents sub-matrix of with the ith data pointed deleted.

To assess the effect of the ith individual on the ML estimate of θ, we consider the following generalised Cook's distance where . Similar to the definition of likelihood displacement (Zhu & Lee, 2001), we define Q-function-based displacement as which quantifies the change of Q-function after deleting the ith observation . Generally, when or are large, the ith observation has a relatively large effect on the ML estimate of θ.

To compute and , we require first evaluating , which indicates that it is necessary to re-run the above presented ECME algorithm for each of n observations. Clearly, it is computationally intensive when sample size n is large. To address the issue, we consider the following one-step approximation of : (35) where is the one-step approximation of , and . Proof of Equation (35(35) ) see Zhu and Lee (2001). Using to replace in and leads to the following approximations and of and : (36) (QD1) (37) (D1) respectively.

To evaluate and , it is necessary to calculate , , and . But all the four quantities have no closed forms for our considered model. To this end, the Monte Carlo approximation is employed. Let be observations of and ω drawn from the joint conditional distribution via the above presented MH algorithm within the Gibbs sampler. Thus, and , and can be approximated by (38) (QQ1) (39) (QQ2) (40) (QQ3) (41) (QQ4) respectively.

The main steps in implementing the above proposed case deletion influence analysis are summarised as follows.

Step (A). Based on observations generated from the joint conditional distribution via the above presented MH algorithm within the Gibbs sampler, calculate , , , and using the approximations given in Equations (35(35) ), (38(38) (QQ1) )–(41(41) (QQ4) ) for.

Step (B). Calculate and via the formulae (36(36) (QD1) ) and (37(37) (D1) ) for .

Step (C). Identify influential observations by inspecting index plots of and based on the corresponding benchmarks, respectively.

5.3. Bayesian local influence analysis

Section 5.1 discusses local influence analysis of EFNSEMs with MNAR based on the Q function of the EM algorithm, but it does not consider perturbation to prior of θ. To this end, we develop a Bayesian local influence approach to assess the effect of simultaneously minor perturbations to the data, the prior and the sampling distribution on the posterior quantities of interest in EFNSEMs with MNAR.

(1) Bayesian perturbation model and manifold

Similar to Zhu et al. (2011), we define a class of perturbation models for simultaneously perturbing to the data, the prior and the sampling distribution as and , where in which , and represent perturbations to the prior, the data and the sampling distribution, respectively. Denote . It is assumed that represents no perturbation.

Following Zhu et al. (2011), the perturbed model can be regarded as a -dimentional manifold under some regularity conditions. The tangent space of at each is composed of functions , where , , and is the kth component of υ. Thus, under some regular conditions, the quantities (42) form the metric tensor of (denoted by , where denotes the expectation taken with respect to the joint probability density function . The th component of measures the amount of perturbation characterised by , while reflects the amount of association between and . If is a diagonal matrix, components of υ are orthogonal to each other and the corresponding perturbation is called an appropriate perturbation. If is not a diagonal matrix, we take the following new perturbation vector such that computed at is equal to , where c is a positive scalar. As an illustration, we consider the following perturbations.

Experiment 1. Consider perturbation to the priors of u, , , , , Φ and ϕ by assuming , , , , , , , and , where , , , , , and are positive scalars, and . Thus, , and indicates no perturbation. The perturbation model constitutes a Riemannian manifold. The tangent space of is consist of Moreover, it is easily shown that

Experiment 2. Consider the perturbation to by assuming . In this case, , and represents no perturbation. The perturbation model forms a manifold. Denote , in which for . The tangent space of is spanned by n quantities: It is easily shown that , where with , and represents the expectation taken with respect to the joint distribution of .

Experiment 3. In this experiment, we consider the perturbation to the sampling distribution by assuming where , is a normalised constant, is a fixed scalar function having zero mean under . This perturbation scheme can be intuitively explained as the perturbation to mean and variance of the normal distribution for , i.e., . In this case, represents no perturbation. The tangent space of is spanned by n quantities: for . It is easily shown that where in which represents the expectation taken with respect to the joint distribution .

Experiment 4. Consider the perturbation to missing data mechanism model (12(12) (LoG) ) by assuming where is defined in Experiment 2, and . Thus, represents no perturbation. The tangent space of the manifold is spanned by . Again, we have which represents the expectation taken with respect to the joint distribution .

(2) Bayesian Local influence measures

Consider the objective function , such as Bayes factor, φ-divergence and posterior mean distance. Let be a geodesic on with . Taking the first-order Taylor expansion of at t=0 leads to , where , , and .

Following Zhu et al. (2011), when , we define the first-order influence (FI) measure in the direction as where and is a positive semi-definite matrix. For the appropriate perturbation , in the unit direction at can be rewritten as Similar to Poon and Poon (1999), we define the first-order adjusted influence measure in the unit direction h at as where with . From the above definition, it is easily seen that to evaluate the first-order adjusted local influence measure , we only require calculating and and selecting an appropriate matrix .

Following the argument given in Section 5.1, for , where is the jth diagonal element of matrix . Again, is used as a benchmark, where and are the mean and standard error of , respectively.

Experiment 5. Consider the logarithm of Bayes factor for comparing υ with on the basis of the observed data : , where . If we take , we have and when . Again, it is rather difficult to calculate . To address the issue, the Monte Carlo approximation method is here adopted. To wit, can be approximated by where observations are generated from the joint posterior distribution via the MH algorithm within the Gibbs sampler.

If , it follows from the Taylor's series expansion that , where and . Following Zhu et al. (2011), The second-order influence measure (SI) in the unit direction is defined as In particular, for an appropriate perturbation , reduces to Thus, the second-order adjusted influence measure at in the unit direction h is defined as where in which . Similarly, diagonal elements of matrix can be used to identify the potential influential observations, misspecified or inappropriate sampling distributions, misspecified missingess data mechanism and misspecified priors. Moreover, is used as a benchmark, where and are the mean and standard error of , respectively.

Experiment 6. We take the objective function as the φ-divergence between two posterior probability density functions before and after introducing perturbation υ, which is defined by where , , and is a convex function with such as the Kullback-Leibler divergence. It is easily shown that and where and is the expectation taken with respect to the joint posterior distribution . In practice, can be approximated by where observations are generated from the joint posterior distribution via the MH algorithm within the Gibbs sampler.

Experiment 7. Define and as the posterior means of before and after introducing υ, respectively, where is some known function of θ. We define Cook's posterior mean distance for characterising the effect of υ on the posterior mean of as where . If we take , we have and , where , which can be approximated by where observations are generated from the joint posterior distribution via the MH algorithm within the Gibbs sampler.

For the above considered perturbation schemes and objective functions, the following four key steps are employed to conduct Bayesian influence analysis.

Step 1. Construct a Bayesian perturbation model .

Step 2. Take an objective function , and calculate , and .

Step 3. If , we calculate the first-order adjusted influence (FIC) measure . However, when , we calculate the second-order adjusted influence (SIC) measure.

Step 4. If (or ), the jth observation is detected as influential, where for k=F and s.

5.4. Bayesian case deletion measures

In this subsection, we consider the effect of deleting a set of observations (denoted as S) on posterior inferences on θ from a Bayesian point of view. In fact, this problem has been investigated in the frequientiest framework, for example, see Section 5.2. A subscript ‘[S]’ represents the relevant quantity with all observations in the set S deleted. Let and be subsamples of and consisting of all the observations in the set S, respectively. The posterior distribution of θ for the observed data set is given by.

Following Zhu et al. (2012), to quantify the effect of deleting all the observations in the set S on the posterior distribution of θ, we define the following Bayesian case influence measure: where , and is a convex function with . is an measure of the distance between two posterior distributions and , and the observations in the data set S are identified as a set of influential observations when is large. Generally, can take various forms, for example, for , for and for corresponding to the Kullback-Leibler divergence (K-L divergence), corresponding to the symmetric K-L divergence, and corresponding to the-distance and the -divergence, respectively.

To measure the effect of deleting all the observations in data set S on the posterior mean of θ, we define the second type of Bayesian case influence measure as (43) where is a user-specified positive definite matrix, and are the posterior means of θ based on the observed data sets and , respectively. Generally, we take as the inverse of the posterior covariance matrix of θ. A large value of indicates that deleting all the observations in data set S leads to a large effect on the posterior mean of θ. Thus, the observations in the data set S are identified as influential.

Although and measure the effect of the data set S on the posterior quantities, their statistical explanations are not the same. assesses the effect of the data set S on the overall posterior distribution, including shape, mode and mean, and is sensitive to posterior distribution. But only measures the effect of the data set S on the posterior mean of θ, and is sensitive to the change of posterior mean.

To compute and , we require calculating the posterior densities and . It is rather difficult to compute them because of latent variables and missing manifest variables involved. To this end, it is desirable to develop an algorithm or a formula computed easily. We define the ratio of the observed data likelihood of θ with and without the data set S as It is easily shown that (44) which indicates that . Thus, can be written as (45) where is the expectation taken with respect to the posterior distribution . In particular, for the K-L divergence (i.e., , is expressed as .

It follows from Equation (44(44) ) and the definitions of and that and , which indicates that can be evaluated by using MCMC samples from the full posterior distribution . In particular, can be computed by averaging the MCMC samples of θ and can be approximated by the inverse of the posterior covariance matrix of the MCMC samples of θ. From the above discussion, it is easily known that computing Equations (43(43) ) and (45(45) ) heavily depend on the computation of.

According to the definition of , we have which shows that it is very cumbersome to approximate . Thus, it is rather unfeasible to compute and because of an intractable multiple integral involved. Thus, some numerical methods such as the trapezoidal rule can be adopted to approximate.

On the other hand, similar to Zhu et al. (2012), one can develop a computationally feasible first-order approximation to and . In particular, under some regularity conditions, can be approximated by and the one-step approximation for is given by where , at . It follows from that which show that and can be approximated by where observations are drawn from the conditional distribution .

6. Model Selection of EFNSEMs with MNAR

6.1. EM-based model selection criterion

By the discussion of Section 4.1, the EM algorithm includes the E-step and the M-step. At the κth iteration of the EM algorithm with the current value of θ, the E-step requires evaluating the Q-function, which is defined as (46) where represents the expectation taken with respect to the conditional distribution , and (47) is called the H-function. The M-step is to maximise to evaluate . When the EM algorithm converges, we can obtain the following three quantities: , , and observations drawn from the joint conditional distribution . According to the above discussion, has no closed expression. Hence, it is impossible to directly evaluate . But, it follows from (46(46) ) that , which shows that can be computed from the Q-function and the H-function when the EM algorithm converges. Thus, inspired by AIC/BIC criterion, we consider the following model selection criterion: (48) where is a penalty function of the data and the fitted model. Different expressions of correspond to different criteria, for instance, when , where d is the dimension of θ, we obtain the AIC; when , the criterion reduces to the BIC. Generally, is selected by users.

It is easily seen from the definition of that is not a direct byproduct of the EM algorithm, and has no closed form. To address the issue, a truncated Hermite approximation is adopted to approximate . The details see Ibrahim et al. (2008).

Since it is rather cumbersome to evaluate the approximation to the integrand of the H-function, it is desirable to present a model selection criterion that does not involve the H-function and only depends on the byproducts of the EM algorithm. To this end, we consider the following model selection criterion by dropping the H-function from (48(48) ): which can be regarded as a reduced form of . When , reduces to the criterion . When , reduces to the criterion .

As Ibrahim et al. (2008) pointed out, the advantage of using is that it is more computationally feasible than , and the disadvantage of is that it results in a criterion with poor model selection properties in some settings where the missing-data proportion is high. Generally, we recommend the usage of over.

6.2. Penalised-EM-based model selection criterion

The above presented EM-based algorithm model selection criterion often suffer from instability and computationally infeasible when the number of the latent variables and covariates is large. To address the issue, various penalisation-based approaches have been proposed to simultaneously estimate parameters and select important latent variables and covariates over the past years. Recently, Tang and Tang (2018) extended the penalised-based approaches to generalised partially nonlinear models with nonignorable missing responses. Here, we extend the penalised-based approaches to EFNSEMs with nonignorable missing data.

Similar to Tang and Tang (2018), at the κth iteration of the EM algorithm with the current value of θ, the E-step requires evaluating the following penalised Q-function: where and are the numbers of different parameters in and , respectively, , and , and are the th, th and th components of Φ, , is the penalty parameter corresponding to the jth element in for , is the penalty parameter corresponding to the jth element in for, is the penalty parameter corresponding to the jth component in Φ for j=1,2,3, while is the penalty parameter corresponding to the jth element in ϕ for , and , , , and , where , , , and .

The M-step is to maximise , which is a quite difficult task in that is a non-differentiable and non-concave function of θ. To address the issue, following Fan and Li (2001), we maximise the second-order Taylor expansions of , , and at , , and , respectively. Thus, the problem of maximising , , and becomes an optimisation problem of the penalised weighted least squares regression, which can be conducted via some appropriate optimisation algorithms such as the local quadratic approximation algorithm (Fan & Li, 2001), local linear approximation algorithm (Zou & Li, 2008). Note that the estimators obtained by maximising with respect to , with respect ton , with respect to Φ and with respect to ϕ, respectively, are not the maximiser of with respect to θ. To address the issue, we here present a hybrid algorithm combining the local linear approximation algorithm and the expectation conditional maximisation algorithm (Meng & Rubin, 1993) to find such that . Iterating the above introduced E-step and the adjusted M-step until the convergence of the EM algorithm leads to the maximum penalised likelihood estimation of θ.

To guarantee that the obtained estimator of θ possesses the well-known oracle property, we require selecting an appropriate penalty parameter λ. Generally, the generalised cross-validation (GCV) and BIC can be employed to select λ. But it is quite difficult to evaluate these criteria in the presence of missing data because of intractable multiple integral involved. Here the above presented criterion is used to choose the penalty parameter λ. To wit, the optimal penalty parameter λ can be evaluated by minimising where is defined in (48(48) ).

Without loss of generality, we assume , where and correspond to the nonzero and zero components of θ, respectively. Let be the true value of θ. Thus, . We write the corresponding decomposition of as , where and correspond to the nonzero and zero components of θ. Similar to Tang and Tang (2018), under some regularity conditions, we can show

1. (Consistency) as ;

2. (Sparsity) Pr(.

6.3. Resampling-based variable selection procedure

Following Long and Johnson (2015), we present a resampling approach combining bootstrap imputation and stability selection to select variables in EFNSEMs with nonignorable missing data.

First, we conduct bootstrap imputation on before variable selection, which is given as follows. We generate bootstrap data sets with the corresponding missing indicator sets based on the observed data set ; and then conduct imputation for each bootstrap data set and using a selected imputation method, which leads to the imputed data set .

Second, given the bootstrap imputed data sets , we consider the following stability selection procedure. We evaluate by maximising for the kth bootstrap imputed data set ; repeat the above procedure for all bootstrap imputed data sets and obtain . Let be the support of (the set of non-zero parameter estimates), which is also called as the estimated active set. The final estimated active set is defined as , where , and is a threshold and is taken as a value between 0.6 and 0.9. Based on the selected active set , we calculate the final parameter estimates via , where is computed by maximising for the kth bootstrap imputed data set , which is the partial data-set of corresponding to .

7. Conclusions and discussions

This paper reviews statistical inference on mean functionals and EFNSEMs in the presence of missing response not at random. Two popular approaches including the EL method, and semiparametric method have been introduced to estimate the mean functionals in the presence of missing response not at random. EM algorithm, Bayesian approach, and Bayesian EL method are presented to estimate parameters in EFNSEMs with nonignorable missing data. Local influence analysis and case-deletion method are developed to assess the effect of minor perturbations on the quantities of interest, and identify the potential influential observations from the frequentist and Bayesian viewpoints, respectively. Also, three model selection approaches, including EM-based model selection criterion, penalised-EM-based model selection criterion, and resampling-based variable selection procedure, are proposed to select the best model among candidate EFNSEMs with nonignorable missing data. Although we mainly review the developments in EFNSEMs with MNAR, it is quite easy to extend these approaches to other models with MNAR.

Recently, Fang and Shao (2016) presented a new model selection criterion, which is referred to as the penalised validation criterion, in the presence of nonignorable nonresponse with unspecified propensity score using the information included in an instrument. A natural idea is to extend the penalised validation criterion of Fang and Shao (2016) to EFNSEMs in the presence of nonignorable missing data with unspecified propensity score. To our best of knowledge, there is little work done on this topic. Also, the aforementioned four model selection methods mainly focus on the setting that the number of explanatory variables is fixed and less than or equal to sample size. Recently, growing-dimensional data analysis receives a huge attention in the absence of missing data. For example, see Fan and Peng (2004) for the nonconcave penalised likelihood in linear regression models; Lam and Fan (2008) for the profile-kernel likelihood in linear regression models; Zou and Zhang (2009) for an adaptive elastic-net procedure in linear regression models; Li, Peng, and Zhu (2011) for a nonconcave penalised M-estimator in linear regression models; Caner and Zhang (2014) for the least squares based adaptive elastic net estimator in generalised method of moments (GMMs); Leng and Tang (2012) for the penalised EL method in estimating equations when likelihood function is unavailable, which is sensitive to model misspecification. To this end, Tang, Yan, and Zhao (2018) presents a penalised exponentially tilted (PET) likelihood for growing dimensional unconditional moment models in the presence of correlation among variables and model misspecification. Hence, extension of the aforementioned model selection techniques to growing-dimensional models with missing data in the presence of model misspecification and outliers is an interesting topic.

Recently, statistical analysis of high/ultrahigh dimensional data, in which the number of candidate explanatory variables increases at an exponential rate of sample size but only a small number of explanatory variables have a relatively large effect on the response, is widely studied. Usually, a two-steps procedure is adopted to analyze the kind of datasets. To wit, the first step is to compress the important explanatory variables to a small number, and the second step is to use the existing penalised variable selection techniques, such as Lasso, SCAD and Adaptive Lasso, to identify the subset of important explanatory variables. To efficiently implement the first step, many feature screening procedures in the absence of missing data have been proposed. For example, see the sure independence screening (SIS) and iterated sure independence screening (ISIS) method (Fan & Lv, 2008) for linear regression models, the sure independent ranking and screening procedure (Zhu, Li, Li, & Zhu, 2011) for multiple-index models, the nonparametric independence screening method (Fan, Feng, & Song, 2011), the robust rank correlation screening procedure (Li, Zhong, & Zhu, 2012) for semiparametric models, the conditional-correlation-based feature screening procedure for varying coefficient models, the GEE-based screening procedure (Xu & Chen, 2014), the score-test-based screening framework (Zhao & Li, 2014), the Pearson correlation sure independence screening procedure based on the partial residual for PLMs Liu, Lee, Zhao (2016). Also, several model-free screening procedures in the absence of missing data have been developed in recent years, for example, see Li et al. (2012), Cui, Li, and Zhong (2015), among others. When responses are subject to missingness, Ding and Wang (2011) presented a Fusion-Refinement procedure for dimension reduction; Lai, Liu, Liu, and Wan (2017) presented a model-free screening procedure combining the IPW method and the Kolmogorov filter method. But, Lai et al.'s (2017) method requires estimating unknown propensity score function, which indicates that it may be subject to the potential misspecification of the respondent probability. Also, Wang and Li (2018) proposed a missing indicator imputation screening procedure for ultrahigh dimensional data in the presence of missing data, which may suffer from the well-known ‘curse of dimensionality’. However, it is a quite challenging and interesting topic to develop a general model-free approach to implement the first step in the presence of nonignorable missing data.

Although we present several approaches to identify the potential influential observations from the frequentist and Bayesian viewpoints, the aforementioned approaches have been developed for regression models with the fixed dimension of explanatory variables. Hence, it is quite interesting to develop some new approaches to detect the potential influential observations in growing-dimensional or high/ultrahigh dimensional models with nonignorable missing data.

Acknowledgements

The authors are grateful to the Editor and two referees for their valuable suggestions that greatly improved the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the grants from the National Natural Science Foundation of China (Grant No.: 11671349) and the Key Projects of the National Natural Science Foundation of China (Grant No.: 11731101).

Notes on contributors

Niansheng Tang

Niansheng Tang is a Professor in Key Lab of Statistical Modeling and Data Analysis of Yunnan Province, Yunnan University.

Yuanyuan Ju

Yuanyuan Ju is a Ph.D. student in Key Lab of Statistical Modeling and Data Analysis of Yunnan Province, Yunnan University.

Appendix. Expressions of and

where , is the kth component of , is the kth row vector of matrix , is the kth row vector of Φ.