
Simulate \(n=5\,000\)  complete values of the binary \(0/1\)  covariate \(x\)  and binary \(0/1\)  outcome \(y\)  from the following 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 \)  parameters are used throughout to make bias comparable across all simulation settings. This sample size is chosen to minimise the issue of small-sample bias associated with the logistic regression [1]};

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, to reflect strong odds ratios (OR) of observing \(x\)  (OR \(= 4.5\)  and \(0.2\) , respectively). For M4, \(\alpha _{y} = 1.5\)  and \(\alpha _{x} = -1.5\)  are chosen as bias in the three MI methods under evaluation is likely to be apparent with these coefficients predicting missingness in \(x\) . For 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\)  appear to work well;

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

Impute missing values in \(x\)  \(M=50\)  times using standard MI and calibrated-\(\delta \)  adjustment MI in turn;

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

