The Jamshidian-Jalal procedure assesses the MCAR assumption by comparing observations belonging to different missingness patterns. Because the test requires complete vectors, it is applied separately to each completed data set. The missingness-pattern groups are always defined using the original incomplete data.

Let a retained completed data set contain $n$ observations, $p$ variables and $g$ missingness patterns. For pattern $r$, let $n_r$ be its frequency, let $Y_{ri}$ denote completed observation $i$ in that pattern, and let

\[ \overline{Y}_r = \frac{1}{n_r}\sum_{i=1}^{n_r}Y_{ri} \]

be the corresponding pattern-specific mean.

Hawkins' transformation (Hawkins, 1981) begins with the pooled within-pattern covariance matrix

\[ \widehat{\Sigma}_{P} = \frac{1}{n-g} \sum_{r=1}^{g} \sum_{i=1}^{n_r} (Y_{ri}-\overline{Y}_r) (Y_{ri}-\overline{Y}_r)^{\mathsf T}. \]

Thus, variation between the pattern-specific means is removed before the common covariance matrix is estimated. Under MCAR and multivariate normality, the covariance matrices of the pattern groups should be homogeneous.

For observation $i$ in pattern $r$, the squared Mahalanobis distance from the pattern mean is

\[ d_{ri} = (Y_{ri}-\overline{Y}_r)^{\mathsf T} \widehat{\Sigma}_{P}^{-1} (Y_{ri}-\overline{Y}_r). \] Define \[ h_{ri}=n_r d_{ri}. \] Hawkins' transformation is then \[ F_{ri} = \frac{(n-g-p)h_{ri}} {p\left\{(n_r-1)(n-g)-h_{ri}\right\}}. \] Under multivariate normality and equality of the covariance matrices across missingness patterns, $F_{ri}$ has an $F$ distribution with $p$ and $n-g-p$ degrees of freedom. Consequently, the upper-tail probability \[ A_{ri} = \Pr\left\{ F_{p,n-g-p}\geq F_{ri} \right\} \]

should follow a uniform distribution on $(0,1)$ within every missingness pattern.

Uniformity is assessed separately in each pattern using a fourth-order Neyman smooth test (Neyman, 1937; Rayner and Best, 1989). Let $\phi_1,\ldots,\phi_4$ denote the first four orthonormal shifted Legendre polynomials. The fourth order means that the first four components of the orthonormal polynomial expansion are retained. These components detect smooth departures from uniformity in four different directions. For pattern $r$, the statistic is

\[ N_r = \frac{1}{n_r} \sum_{\ell=1}^{4} \left\{ \sum_{i=1}^{n_r} \phi_{\ell}(A_{ri}) \right\}^{2}. \]

For a pattern containing at least $\mathtt{usechisq}$ observations, $N_r$ is compared with a chi-square distribution having four degrees of freedom:

\[ N_r \mathrel{\dot{\sim}} \chi^2_4. \]

For smaller patterns, the reference distribution is approximated by generating nsimul samples of size $n_r$ from the uniform distribution and recomputing the Neyman statistic for each sample.

Let $p_r$ be the Neyman smooth-test $p$-value obtained for pattern $r$.

The pattern-specific results are combined using Fisher's method (Fisher, 1932), giving the statistic

\[ H = -2\sum_{r=1}^{g}\log(p_r). \] Under the joint null hypothesis, \[ H \mathrel{\dot{\sim}} \chi^2_{2g}. \]

A small combined Hawkins $p$-value indicates that at least one pattern departs from the expected uniform behaviour. Such a rejection can be caused either by covariance heterogeneity, which provides evidence against MCAR under the assumptions of the procedure, or by departure from multivariate normality.

When method='hawkins', only this combined Hawkins test is reported. When method='auto', Hawkins' test is computed first. If at least one imputation-specific Hawkins p-value is smaller than alpha, the function also applies the k-sample Anderson-Darling test. When method='nonparametric', the Anderson-Darling test is always computed.

The k-sample Anderson-Darling test of Scholz and Stephens (1987) compares the distributions of the transformed values $F_{ri}$ across the $g$ retained missingness patterns. It extends the weighted empirical-distribution criterion introduced by Anderson and Darling (1952) to the k-sample problem. Its null hypothesis is

\[ H_0:\mathcal{F}_1=\cdots=\mathcal{F}_g, \]

where $\mathcal{F}_r$ denotes the distribution of $F_{ri}$ in pattern $r$. Under MCAR, these pattern-specific distributions should be equal.

Let

\[ z_1<\cdots<z_L \]

denote the distinct ordered values in the pooled sample of transformed observations, excluding the largest pooled value. For each $z_l$, let $h_l$ be its multiplicity in the pooled sample, let $H_l$ be the number of pooled observations not greater than $z_l$, and let $M_{rl}$ be the number of observations from pattern $r$ not greater than $z_l$.

The contribution of pattern $r$ to the Anderson-Darling statistic is

\[ A_r = \frac{1}{n_r} \sum_{l=1}^{L} h_l \frac{ \left(nM_{rl}-n_rH_l\right)^2 }{ H_l(n-H_l) }. \] The complete k-sample statistic is \[ A^2_{g,n} = \frac{1}{n} \sum_{r=1}^{g}A_r. \]

The denominator $H_l(n-H_l)$ gives relatively high weight to distributional differences occurring in the tails. The multiplicities $h_l$ provide an adjustment for tied transformed values.

Under the null hypothesis, the expected value of $A^2_{g,n}$ is approximately $g-1$. The function computes its finite-sample variance, denoted by $\sigma^2_{g,n}$, and forms the standardized statistic

\[ T_{g,n} = \frac{ A^2_{g,n}-(g-1) }{ \sqrt{\sigma^2_{g,n}} }. \]

Large positive values of $T_{g,n}$ indicate stronger differences among the pattern-specific distributions and therefore correspond to small Anderson-Darling $p$-values.

Following the approximation of Scholz and Stephens (1987), the Anderson-Darling $p$-value is not computed from a chi-square or an $F$ reference distribution. The implementation uses critical values corresponding to the five upper-tail probabilities

\[ p_l \in \{0.25,0.10,0.05,0.025,0.01\}. \]

For each probability $p_l$, the corresponding critical value is approximated as a function of the number of retained patterns:

\[ q_l(g) = b_{0l} + \frac{b_{1l}}{\sqrt{g-1}} + \frac{b_{2l}}{g-1}. \]

The tabulated probabilities are transformed to log-odds:

\[ c_l = \log\left( \frac{1-p_l}{p_l} \right). \] A cubic spline interpolates the value $c$ corresponding to the observed standardized statistic $T_{g,n}$. The approximate upper-tail probability is then \[ p_{\mathrm{AD}} = \frac{1}{1+\exp(c)}. \]

Therefore, $\mathtt{out.ADpvalue}(m)$ is the approximate Anderson-Darling p-value obtained from completed data set $m$. A small value indicates that the distributions of the Hawkins-transformed quantities differ across the retained missingness patterns.

Because the calculation uses spline interpolation, and may use spline extrapolation when $T_{g,n}$ lies outside the tabulated range, very small or very large Anderson-Darling p-values should be interpreted as approximate tail probabilities.

The k-sample Anderson-Darling test is distribution free and is used to help distinguish a Hawkins rejection caused mainly by nonnormality from a rejection caused by differences among the missingness-pattern distributions. In the automatic procedure, a significant Anderson-Darling result provides evidence against MCAR within the MCAR-versus-MAR framework. A nonsignificant Anderson-Darling result following a significant Hawkins test suggests that nonnormality is the more plausible explanation.

When $M$ completed data sets are analysed, the entire procedure is repeated separately for each imputation. The imputation-specific statistics and p-values are retained in $\mathtt{out.hawkinsStat}$, $\mathtt{out.hawkinsP}$, $\mathtt{out.ADstat}$ and $\mathtt{out.ADpvalue}$.

The main reported statistic and $p$-value are the medians of the imputation-specific quantities associated with the selected branch:

\[ T_{\mathrm{reported}} = \mathop{\mathrm{median}}_{m=1,\ldots, M} T^{(m)}, \qquad p_{\mathrm{reported}} = \mathop{\mathrm{median}}_{m=1,\ldots, M} p^{(m)}. \] In particular, when the Anderson-Darling branch supplies the main result, \[ \mathtt{out.pvalue} = \mathop{\mathrm{median}} \left\{ p_{\mathrm{AD}}^{(1)},\ldots, p_{\mathrm{AD}}^{(M)} \right\}. \]

This median is a descriptive summary across imputations. It is not a Rubin-rules or Meng-Rubin pooled p-value. The complete vectors of imputation-specific p-values should therefore also be examined.

The global Hawkins and Anderson-Darling statistics can also be decomposed into nonnegative pattern-specific contributions. For Hawkins' test, let $p_r^{(m)}$ be the Neyman smooth-test p-value for retained pattern r in completed data set m. Its contribution is

\[ C_{\mathrm{H},r}^{(m)} = -2\log\left(p_r^{(m)}\right), \] and therefore \[ H^{(m)} = \sum_{r=1}^{g}C_{\mathrm{H},r}^{(m)}. \] The mean Hawkins contribution returned for pattern r is \[ \overline{C}_{\mathrm{H},r} = \frac{1}{M} \sum_{m=1}^{M} \left[-2\log\left(p_r^{(m)}\right)\right]. \] For the Anderson-Darling test, let $C_{\mathrm{AD},r}^{(m)}$ denote the normalized group contribution returned for pattern r in completed data set m. These contributions satisfy \[ A_{g,n}^{2(m)} = \sum_{r=1}^{g}C_{\mathrm{AD},r}^{(m)}, \] and their mean is \[ \overline{C}_{\mathrm{AD},r} = \frac{1}{M} \sum_{m=1}^{M}C_{\mathrm{AD},r}^{(m)}. \] Consequently, the sums of the returned mean contributions satisfy \[ \sum_{r=1}^{g}\overline{C}_{\mathrm{H},r} = \frac{1}{M}\sum_{m=1}^{M}H^{(m)}, \qquad \sum_{r=1}^{g}\overline{C}_{\mathrm{AD},r} = \frac{1}{M}\sum_{m=1}^{M}A_{g,n}^{2(m)}. \]

Thus, the mean contributions decompose the mean imputation-specific statistic. They do not generally decompose $\mathtt{out.stat}$, because the latter is the median of the imputation-specific statistics.

To construct the Pareto plot, let

\[ \overline{C}_{(1)} \geq \overline{C}_{(2)} \geq \cdots \geq \overline{C}_{(g)} \] denote the mean contributions sorted from largest to smallest. The cumulative percentage after the first r sorted patterns is \[ 100 \frac{\sum_{j=1}^{r}\overline{C}_{(j)}} {\sum_{j=1}^{g}\overline{C}_{(j)}}. \]

Large bars identify patterns that contribute most strongly, on average across completed data sets, to the corresponding global statistic. A large contribution is a diagnostic indication and should not be interpreted as a separate formal rejection for that pattern.

Option nsimul affects only the Monte Carlo calibration of the pattern-specific Neyman tests in the Hawkins branch. It does not affect the Anderson-Darling p-values.

When imputed is empty, mdJJtest estimates location and covariance using mdEM and generates stochastic completed data sets using mdImputeStochastic. This is a joint multivariate-normal imputation and is not identical to the default chained normal-regression imputation used by mice::mcar. To compare the MATLAB and R implementations using exactly the same completed data sets, supply them through option imputed.

The procedure addresses the MCAR-versus-MAR framework. It cannot distinguish MCAR from MNAR. Consequently, a nonsignificant result does not establish MCAR and does not rule out an MNAR mechanism.