
For cases 2 and 3, the following two steps are performed to incorporate the sampling behaviour of \(\hat{p}_{x}^{\text{pop}}\) , which is estimated in an external dataset of size \(n^{\text{ex}}\) , into the data generating mechanism in repeated simulations.

Simulate \(n^{\text{ex}} = 10\,000\)  (case 2) or \(1\,000\)  (case 3) complete values of the binary \(0/1\)  covariate \(x\)  from the model
\(x \sim \text{Bernoulli}\left(p_{x}^{\text{pop}} = 0.7\right);\) 

Obtain the sample proportion \(\hat{p}_{x}^{\text{pop}}\)  of \(x\) , which is an unbiased estimate of the population proportion \(p_{x}^{\text{pop}}\) ;


Simulate \(n=5\,000\)  complete values of the binary \(0/1\)  covariate \(x\)  and binary \(0/1\)  covariate \(y\)  from the models
\(&x \sim \text{Bernoulli}\left(p_{x}^{\text{pop}} = 0.7\right); \nonumber \\&\text{logit}\left[p\left(y=1 \mid x\right)\right] = \beta _{0} + \beta _{x}x ,\) 

where \(\beta _{0}\)  and \(\beta _{x}\)  are arbitrarily set to \(\text{ln}\left(0.5\right)\)  and \(\text{ln}\left(1.5\right)\) , respectively. The same values of the \(\beta \)  coefficients are used throughout to make bias comparable across all simulation settings;

Simulate a binary indicator of response \(r\)  of \(x\)  from each of the selection models M1–M4 (Table REF ). Values of \(1.5\)  and \(-1.5\)  are chosen for \(\alpha _{y}\)  and \(\alpha _{x}\)  in M2 and M3, respectively. For M4, \(\alpha _{y} = 1.5\)  and \(\alpha _{x} = -1.5\)  are used. In all selection models, \(\alpha _{0}\)  is altered to achieve approximately \(45\%\)  missing \(x\) . For M1, \(\alpha _{0}\)  is calculated directly as \(\text{ln}\left(\frac{0.55}{0.45}\right)\) ; for M2–M4, \(\alpha _{0} = -0.2; 1.35\) ; and \(0.75\)  are used;

For \(i=1, \ldots , 5\,000\) , set \(x_{i}\)  to missing if \(r_{i} = 0\) ;

Impute missing values in \(x\)  \(M=10\)  times using standard MI and calibrated-\(\delta \)  adjustment MI in turn. For cases 2 and 3, calibrated-\(\delta \)  adjustment MI is performed as follows.

Draw a value \(\tilde{p}_{x}^{\text{pop}}\)  from the normal approximation \(N\left(\hat{p}_{x}^{\text{pop}}, \frac{\hat{p}_{x}^{\text{pop}}\left(1-\hat{p}_{x}^{\text{pop}}\right)}{n^{\text{ex}}}\right)\) , with values of \(n^{\text{ex}} = 10\,000\)  (case 2) and \(1\,000\)  (case 3). This is done by first taking a draw \(\tilde{z}\)  from the standard normal distribution, \(z \sim N\left(0,1\right)\) , followed by drawing \(\tilde{p}_{x}^{\text{pop}} = \hat{p}_{x}^{\text{pop}} + \tilde{z}\sqrt{\frac{\hat{p}_{x}^{\text{pop}}\left(\hat{p}_{x}^{\text{pop}}\right)}{n^{\text{ex}}}}\) ;

Derive the calibrated-\(\delta \)  adjustment and perform MI according to the algorithm set out in subsubsec3.1.1, using \(\tilde{p}_{x}^{\text{pop}}\)  as the reference proportion;


For each MI method, fit the analysis model (REF ) to each completed dataset and combine the results using Rubin's rules.[1]}, [2]}

