This function implements the three diagnostics from Bojinov, Pillai and Rubin:

$\mathbf{ccm}$: comparison of conditional means using nested Gaussian linear models and Bonferroni correction.

$\mathbf{dtmm}$: direct testing of the missingness mechanism through logistic models and multiple imputation. The likelihood-ratio statistics are combined using the Meng-Rubin $D_L$ procedure, often referred to as the $D_3$ procedure in the multiple-imputation literature.

​ $\mathbf{cop}$: gaussian-copula diagnostic based on posterior conditional covariances between partially observed variables and missingness indicators.

---------------------------------------------------.

$\mathbf{ccm}$ method.

Let $J_m$ denote the set of variables that contain missing values and $J_f$ the set of fully observed variables, and let

\[ q=|J_m| \]

denote the number of partially observed variables. Under the assumptions considered by Bojinov, Pillai and Rubin, a missingness indicator $R_k$ should not depend on another partially observed variable $Y_j$ after conditioning on $Y_{J_f}$.

The comparison-of-conditional-means procedure tests this implication by comparing, for every $j$ different from $k$, the nested linear models

\[ Y_j \sim Y_{J_f} \] and \[ Y_j \sim Y_{J_f} + R_k + Y_{J_f}:R_k. \]

$\mathbf{dtmm}$ method.

The direct procedure (dtmm) tests the postulated missingness mechanism separately for each missingness indicator $R_k$. For a given $R_k$, the reduced logistic model contains only an intercept and the fully observed variables $Y_{J_f}$:

\[ \mathrm{logit}\{ \Pr(R_k=1\mid Y_{J_f}) \} = \alpha_k+Y_{J_f}^{\mathsf T}\beta_k. \] The full model additionally contains all the partially observed variables after imputation. Thus, in completed data set $m$, \[ \mathrm{logit}\{ \Pr(R_k=1\mid Y_{J_f},Y_{J_m}^{(m)}) \} = \alpha_k^{(m)} +Y_{J_f}^{\mathsf T}\beta_k^{(m)} +(Y_{J_m}^{(m)})^{\mathsf T}\gamma_k^{(m)}. \]

The response $R_k$ and the fully observed variables do not change across imputations. Consequently, the reduced model is identical in all completed data sets and is fitted only once. The full model, on the other hand, depends on the imputed values and is fitted separately in each of the $M$ completed data sets.

Let $\widehat\theta_k^{(m)}$ denote the vector of full-model coefficient estimates obtained from completed data set $m$. The Meng-Rubin $D_L$ procedure first forms the common pooled coefficient vector

\[ \overline\theta_k = \frac{1}{M} \sum_{m=1}^M \widehat\theta_k^{(m)}. \] Each completed data set is then reevaluated at this same pooled coefficient vector, without refitting the model. If $D_{F,m}(\widehat\theta_k^{(m)})$ is the full-model deviance evaluated at its own maximum-likelihood estimate and $D_{F,m}(\overline\theta_k)$ is the deviance evaluated at the pooled coefficients, their average difference \[ \frac{1}{M} \sum_{m=1}^M \left\{ D_{F,m}(\overline\theta_k) - D_{F,m}(\widehat\theta_k^{(m)}) \right\} \]

measures the loss of fit caused by forcing all imputations to use a common full-model parameter vector. It therefore measures the disagreement among the imputation-specific estimates and is used to estimate the additional uncertainty due to missing information.

Since the full model adds the $q$ partially observed variables, the estimated relative increase in variance is

\[ r_L = \frac{M+1}{q(M-1)} \frac{1}{M} \sum_{m=1}^M \left\{ D_{F,m}(\overline\theta_k) - D_{F,m}(\widehat\theta_k^{(m)}) \right\}. \] Let $D_R$ denote the deviance of the reduced model. The numerator used by this function is \[ d_L = \max\left\{0, D_R - \frac{1}{M} \sum_{m=1}^M D_{F,m}(\overline\theta_k) \right\}. \] The maximum with zero is only a numerical safeguard against tiny negative values caused by finite-precision calculations. The pooled Meng-Rubin statistic used by this function is then \[ D_L = \frac{d_L}{q(1+r_L)}. \]

Thus $d_L$ measures the improvement of the full model over the reduced model after imposing common pooled coefficients, whereas $r_L$ quantifies the additional uncertainty arising when the imputation-specific full-model estimates disagree substantially.

The relative increase $r_L$ can also be expressed through the corresponding fraction of missing information,

\[ \lambda_L=\frac{r_L}{1+r_L}. \]

This transformation maps $r_L\in[0,\infty)$ into $\lambda_L\in[0,1)$. Values close to zero correspond to little missing-information uncertainty, whereas values approaching one correspond to a large missing-information component. The DTMM diagnostic plot therefore uses $\lambda_L$ on the vertical axis, while the original $r_L$ is reported in out.DLtable.

The resulting statistic is compared with the finite-$M$ $F$ reference distribution of Meng and Rubin (1992).

A small $p$-value indicates that, after conditioning on the fully observed variables, at least one partially observed variable provides additional information about the missingness indicator $R_k$.

$\mathbf{cop}$ method.

The Gaussian-copula procedure uses the extended rank-likelihood sampler of Hoff (2007). For every retained posterior correlation matrix, the function computes the covariance between the partially observed outcomes and their missingness indicators conditional on the fully observed variables. A component is rejected when its Bonferroni-adjusted credible interval does not contain zero.

The original R implementation uses mice for method 'dtmm'. When option imputed is empty, this MATLAB implementation instead uses the joint-normal EM and stochastic-imputation functions already available in FSDA.

$\large{\text{PLOTS IN OUTPUT}}$.

If plots is true, the graphical summary depends on the selected diagnostic method:

$\mathbf{ccm}$ Produces a q-by-q heatmap of the pairwise F-test p-values. Rows correspond to missingness indicators $R_i$ and columns to partially observed variables $Y_j$. Diagonal cells are not tested. Small p-values are shown using reddish colours, and cells significant after the Bonferroni correction are displayed in red.

$\mathbf{dtmm}$ Produces a bubble scatter plot that separates the two components entering the Meng-Rubin $D_L$ statistic,

\[ D_L=\frac{d_L}{q(1+r_L)}. \]

The horizontal coordinate $d_L/q$ measures the likelihood-ratio improvement in fit per tested restriction obtained by augmenting the reduced missingness model with the partially observed variables.

More precisely, $d_L$ is the reduction in deviance obtained by moving from the reduced model to the full model when the latter is evaluated using the common coefficient vector pooled across imputations.

The vertical coordinate is

\[ \lambda_L=\frac{r_L}{1+r_L}, \]

a bounded representation of the missing-information component. There is one bubble for each missingness indicator, labelled with the corresponding variable name. Bubble size and colour increase as the associated $p$-value decreases; a black ring identifies components rejected at the Bonferroni-adjusted significance level.

The vertical reference line is the Bonferroni-adjusted threshold for $d_L/q$ under $r_L=0$. The horizontal reference line $\lambda_L=0.5$, equivalently $r_L=1$, is an interpretative benchmark for the missing-information component. These reference lines are graphical guides; the formal decision is based on the finite-$M$ Meng-Rubin $p$-value.

For graphical readability, numerator values exceeding six times the reference threshold are displayed at that upper bound and marked by a right-pointing symbol. Their exact values remain available in out.DLtable.

$\mathbf{cop}$ Produces a q-by-q heatmap of the two-sided posterior tail probabilities for the latent conditional covariances between missingness indicators $R_i$ and partially observed variables $Y_j$, conditional on the fully observed variables. Rows correspond to $R_i$ and columns to $Y_j$. Diagonal cells are not tested. Small posterior tail probabilities are shown using reddish colours, and rejected components are displayed in red.

The plot is created in a new MATLAB figure. Setting plots to false suppresses all graphical output but does not affect the numerical results returned in out.