Longitudinal data play an important role in many biomedical and social problems. Multivariate longitudinal data is a generalization of univariate longitudinal data where there are many dependent outcomes obtained at many occasions. Missing data frequently occur in longitudinal studies for many reasons such as subjects’ moving or medical issues. Missingness could affect the bias and precision of parameter estimations. In this paper, we fill this gap with a dedicated simulation to study the missingness effects on the parameter estimations for multivariate longitudinal data. The simulation study is conducted to evaluate the effects of the correlation in the parameter estimations for missing data mechanisms, that are MCAR, MAR and MNAR. Also, several analysis methods for handling missing data are used to reduce the bias of the parameter estimations when there is missingness in the responses, namely Completed cases (CC), Mean substitution (MS), Last Observation Carried Forward (LOCF) and Regression Imputation (RI). The inverse propensity weighed (IPW) GEE is used for the MAR mechanisms to investigate its performance with dropout missingness. Based on evaluation measurements such as root mean square error (RMSE) and coverage probability (CP), we found our results in agreement with Little¹ who indicated MCAR is the appropriate mechanism for the longitudinal data. We found that severe correlation over the occasions may affect the parameter estimations for both complete and incomplete data. Also, we found the MS handling method, after extension to accommodate the multivariate structure, largely has good results even when one outcome is MAR and the second is MCAR. The inverse propensity weighed GEE shows some good results to treat the dropout in MAR mechanism especially when the correlation is induced over the outcomes.

The nature of longitudinal data is repeated measurements over many occasions. Missingness frequently occurs during a longitudinal study because of circumstances such as a subject moving, medical illness or administrative reasons. In this paper, our focus is multivariate longitudinal data, for which there are more than one outcome measured at many occasions. The missingness in multivariate longitudinal data can happen such that all outcomes are missing at a particular time or partially such that only a subset of outcomes is not observed at that time. Missingness is considered a big issue in statistical analysis since its effects are to reduce the statistical power and increase bias. The statistical analysis should pay attention to the problem of missing data and use some statistical missing data handling methods to reduce the bias and reach better estimates. Missing data mechanisms were classified by Little¹ into three types:

1. Missing completely at random (MCAR), when the missingness is independent of the observed and unobserved responses;
2. Missing at random (MAR), when the missingness may depend on observed responses, but is independent of the unobserved responses;
3. Missing not at random (MNAR), when missingness depends on both observed and unobserved responses. MNAR is often called “informative missingness.” In real life, the missingness mechanism is unknown and it is not a trivial issue to assume that missingness is MCAR or MAR. It is important to assume the appropriate mechanism based on the nature of the study to avoid biased estimates. In this paper, we generated artificial data through a simulation study to perform statistical analysis of incomplete binary multivariate longitudinal data in order to answer many questions. First, are the estimated correlation parameters different for the completed data and incomplete data, controlling the correlation over the outcomes and occasions? Second, are the estimated regression parameters valid for the three cases of missingness in all outcomes, in just partial outcomes, or mixed missing mechanisms for the outcomes. The third question is about the best method for handling missing data among CC, MS, LOCF and RI to treat the incomplete binary multivariate longitudinal data.

Many researchers have conducted simulation studies to investigate the effects of missing data handling methods for univariate longitudinal data. For example Myers² conducted a simulation for longitudinal data to compare the complete case method and multiple imputation methods in clinical trials and found some limitations of using multiple imputation. Also, Touloumi et al.³ designed a simulation study to investigate the impact of missing data due to drop out in longitudinal studies for six methods such as GEE, weighted and unweighted ordinary least square. They found the MCAR assumption is doing well for all their methods while MAR and MNAR generate biased estimates. Newman⁴ compared six missing data techniques based on simulation study. The comparison was for the three mechanisms of the missingness (MCAR, MAR, and MNAR) and for three levels of missingness 25%, 50% and 75%. Their results support a multiple imputation approach. Also, Hening⁵ performed a comparison for five imputation methods for the missing data (mean substitution, median substitution, zero values, hot deck and MI) with 20% missing rates for the first year students retention data in Ohio university. His dissertation indicates that mean imputation and median imputation yield good performance in precision. Ali et al.⁶ investigated four imputation methods (complete case analysis, mean substitution, and multiple imputation with and without inclusion of the outcome in the imputation model) under MCAR, MAR and MNAR based on simulation study to figure out the best approach. They found, based on their model, that the estimates for multiple imputation were least biased and most accurate.

The purpose of this paper is an investigation of the effect of missing data and the performance of missing data handling methods, especially for binary multivariate longitudinal data through a simulation study. Multivariate binary longitudinal data are generated for specified correlation structures and for different missing data mechanisms. The organization of this paper is as follows: section 3 contains the simulation design. It contains the model details, correlation scenarios, missing data mechanisms, handling missing data methods and simulation evaluation measurements. In section 4, the simulation results are produced. Finally, discussion and concluding remarks are in section 5.

The multivariate longitudinal data structure is an extension of the univariate longitudinal structure to more than one outcome. Each individual $i$ has a vector of responses for different outcomes, $k\mathrm{<mo>=</mo>1,2,}\dots \mathrm{,}K$. Also, each individual is measured at different occasions, $j\mathrm{<mo>=</mo>1,2}\dots \mathrm{,}{J}_i$, and has cluster size $n_i\mathrm{<mo>=</mo>}{J}_{i}K$. To simplify these notations, we will refer to $J_i$ as $J$, which is the number of occasions for all individuals. The vector of completed responses for subject $i$ is:

$Y_i\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}{Y}_{i11}\mathrm{,}{Y}_{i21}\mathrm{,}\dots \mathrm{,}{Y}_{iJ1}\mathrm{,}{Y}_{i12}\mathrm{,}{Y}_{i22}\mathrm{,}{Y}_{i32}\mathrm{,}\dots \mathrm{,}{Y}_{iJ2}\mathrm{,............,}{Y}_{i1K}\mathrm{,}{Y}_{i2K}\mathrm{,}\dots \mathrm{,}{Y}_{iJK}{\mathrm{<mo\; stretchy="false">]</mo>}}^T.$

We conducted a simulation study in order to explore the changes of the parameter estimations for different missing data mechanisms when the data structure is multivariate longitudinal binary data. One of the goals to design the simulation study was to control the missing data pattern and the correlation for multivariate structure through the regression model.

Simulation model

We will use generalized estimating equations (GEE), an extension of generalized linear modeling to longitudinal data, to specify the simulation model and for the analysis. Thus, we specify a marginal model for the correlated responses. Estimation via GEE yields consistent regression parameter estimations despite the lack of full likelihood specification when the data are complete. This consistency also holds when data are MCAR, as shown in Little.¹ We fit the GEE model of Shelton et al.⁷ to estimate the effects for each outcomes separately using a Kronecker product approach to account for the correlation structure. Let $X_i\mathrm{<mo>=</mo>0,1}$ be the treatment assignment covariate. The time covariate is $t_j\mathrm{<mo>=</mo>1,2,3}$ for three occasions. Let $Y_{ijk}$ be the binary response that is measured at time $j\mathrm{,}j\mathrm{<mo>=</mo>1,2,3}$ for observation $i\mathrm{,}i\mathrm{<mo>=</mo>1,2,3,}\dots \mathrm{,}N$ and for outcome $k\mathrm{<mo>=</mo>1,2}$. Then, we assume the logistic model

$logit\mathrm{<mo\; stretchy="false">(</mo>}E\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{ijk}\mathrm{<mo>|</mo>}{X}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo>=</mo>}{\beta }_{0k}+{\beta }_{1k}{X}_i+{\beta }_{2k}{t}_j$

with the true parameter values ${\beta }_{01}\mathrm{<mo>=</mo>0.2,}{\beta }_{02}\mathrm{<mo>=</mo>0.1,}{\beta }_{11}\mathrm{<mo>=</mo>0.05,}{\beta }_{12}\mathrm{<mo>=</mo>}-\mathrm{0.15,}{\beta }_{21}\mathrm{<mo>=</mo>0.05,}$ and ${\beta }_{22}\mathrm{<mo>=</mo>}-0.25$. Given the vector of treatment covariate $X_i$, the $Y_{ijk}$ for subject $i$ is assumed to follow the Bernoulli distribution, $Y_{ijk}\mathrm{<mo>|</mo>}{X}_i\sim Ber\mathrm{<mo\; stretchy="false">(</mo>}E\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{ijk}\mathrm{<mo>|</mo>}{X}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$. We did the simulation for the marginal regression model based on specified correlation structures and for different missing data mechanisms. We applied the method of Alzahrani⁸ that used the bridge distribution for the random effect in the mixed model to generate the correlated binary data. We generated $N\mathrm{<mo>=</mo>250}$ samples of correlated binary data, conducted the simulation study for eight missing data mechanisms and for five scenarios (correlation structures) to explore the properties of the model when the correlation is induced over the outcomes and the occasions. In the following subsections, there are more explanations about correlation scenarios, the missing data mechanisms, handling missing data methods and evaluation measurements.

Correlation scenarios

Multivariate longitudinal data has a complicated correlation structure in which the correlation has two factors, over the occasions and over the outcomes. This simulation study explores the effects of increasing correlation over each of these two factors. We start by explaining the nature of the correlation in multivariate longitudinal data. The correlation matrix $R\mathrm{<mo\; stretchy="false">(</mo>}\gamma \mathrm{<mo\; stretchy="false">)</mo>}$ is a function of ${\gamma }_{jk\mathrm{,}{j}^{\mathrm{<mo>\text{'}</mo>}}{k}^{\mathrm{<mo>\text{'}</mo>}}}$, where $\mathrm{<mo>\{</mo>}{\gamma }_{jk\mathrm{,}{j}^{\mathrm{<mo>\text{'}</mo>}}{k}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{<mo>\}</mo>}$ represents the collection of within subject correlation parameters of size $\left(\frac{JK}{2}\right)$ of all non-redundant pairwise correlation parameters. For our setup with $J\mathrm{<mo>=</mo>3}$ and $K\mathrm{<mo>=</mo>2}$, there are 15 correlation parameters given by:

${\gamma }_{jk\mathrm{,}{j}^{\mathrm{<mo>\text{'}</mo>}}{k}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{<mo>=</mo>}\frac{P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{jk}\mathrm{<mo>=</mo>1,}{Y}_{j^{\mathrm{<mo>\text{'}</mo>}}{k}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo>}-P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{jk}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{j^{\mathrm{<mo>\text{'}</mo>}}{k}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo>}}{\sqrt{P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{jk}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{j^{\mathrm{<mo>\text{'}</mo>}}{k}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{jk}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo><mo\; stretchy="false">(</mo>1}-P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{j^{\mathrm{<mo>\text{'}</mo>}}{k}^{\mathrm{<mo>\text{'}</mo>}}}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}}$

$\begin{array}{l}R\mathrm{<mo>=</mo>}{Y}_{11}{Y}_{21}{Y}_{31}{Y}_{12}{Y}_{22}{Y}_{32}{Y}_{11}1{\gamma }_{\mathrm{11,21}}{\gamma }_{\mathrm{11,31}}{\gamma }_{\mathrm{11,12}}{\gamma }_{\mathrm{11,22}}{\gamma }_{\mathrm{11,32}}{Y}_{21}-1\\ {\gamma }_{\mathrm{21,31}}{\gamma }_{\mathrm{21,12}}{\gamma }_{\mathrm{21,22}}{\gamma }_{\mathrm{21,32}}{Y}_{31}-1{\gamma }_{\mathrm{31,12}}{\gamma }_{\mathrm{31,22}}{\gamma }_{\mathrm{31,32}}{Y}_{12}---1{\gamma }_{\mathrm{12,22}}{\gamma }_{\mathrm{12,32}}{Y}_{22}----1{\gamma }_{\mathrm{22,32}}{Y}_{32}-----1\end{array}$

To reduce size of the unknown correlation parameter vector $\gamma$, we assume there are three components that build up the correlation structure $R$ in the multivariate longitudinal data. The inter-outcome ($\upsilon$), the intra-outcome ($\alpha$); and the cross association ($\tau$). The $\alpha s$ represent the correlation between the outcomes in the same occasion. The $\upsilon s$ represent the correlation within outcomes at different occasions, and $\tau s$ represent the correlation between different outcomes measured at different occasions. Thus, $R$ can be written as follows:

$R\mathrm{<mo>=</mo>}{Y}_{11}{Y}_{12}{Y}_{13}{Y}_{21}{Y}_{22}{Y}_{23}{Y}_{11}1{\upsilon }_1{\upsilon }_2{\alpha }_1{\tau }_1{\tau }_2{Y}_{12}-1{\upsilon }_3{\tau }_3{\alpha }_2{\tau }_4{Y}_{13}--1{\tau }_5{\tau }_6{\alpha }_3{Y}_{21}---1{\upsilon }_4{\upsilon }_5{Y}_{22}----1{\upsilon }_6{Y}_{23}-----1$

We conducted five scenarios under the assumption of exchangeability for each correlation type $\upsilon \mathrm{,}\alpha$ and $\tau$ and $\tau \mathrm{<mo>=</mo>0}$ for all scenarios. We faced a difficulty to do the rest of the scenarios to the restrictions on generating correlated binary outcomes from the regression model for specified correlation strictures. The following table shows the values for each correlation parameter in each scenario for the correlation matrix $R$ (Table 1).

Missing data mechanisms

It is important to investigate the nature of the missing data to get valid inference. The missing data mechanism can be defined as the probability distribution of the missing indicator variable $R_i\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>0,1<mo\; stretchy="false">)</mo>}$ that takes value 1 when the response is missing, and 0 if not. The vector of complete responses for subject $i$ is:

$Y_i\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}{Y}_{i11}\mathrm{,}{Y}_{i21}\mathrm{,}\dots \mathrm{,}{Y}_{iJ1}\mathrm{,}{Y}_{i12}\mathrm{,}{Y}_{i22}\mathrm{,}{Y}_{i32}\mathrm{,}\dots \mathrm{,}{Y}_{iJ2}\mathrm{,............,}{Y}_{i1K}\mathrm{,}{Y}_{i2K}\mathrm{,}\dots \mathrm{,}{Y}_{iJK}{\mathrm{<mo\; stretchy="false">]</mo>}}^T$

Let $R_i$ be the vector of response missingness indicators,

$R_i\mathrm{<mo>=</mo><mo\; stretchy="false">[</mo>}{R}_{i11}\mathrm{,}{R}_{i21}\mathrm{,}\dots \mathrm{,}{R}_{iJ1}\mathrm{,}{R}_{i12}\mathrm{,}{R}_{i22}\mathrm{,}{R}_{i32}\mathrm{,}\dots \mathrm{,}{R}_{iJ2}\mathrm{,............,}{R}_{i1K}\mathrm{,}{R}_{i2K}\mathrm{,}\dots \mathrm{,}{R}_{iJK}{\mathrm{<mo\; stretchy="false">]</mo>}}^T\mathrm{,}$

with $R_{ijk}\mathrm{<mo>=</mo>1}$ when $Y_{ijk}$ is not observed, and $R_{ijk}\mathrm{<mo>=</mo>0}$ when $Y_{ijk}$ is observed. In this paper, we do not consider missingness in the covariates. Given $R_i$, the complete set of responses $Y_i$ can be divided into two groups: $Y_i^O$ and $Y_i^M$. $Y_i^O$ denotes the vector of observed responses and $Y_i^M$ the vector of missing responses. Setting up the three mechanisms is done by the approaches in Gibbons.⁹ Assuming the responses at the first time are fully observed, then the three missing mechanisms are as follows.

Missing completely at random MCAR:

The missing data pattern is considered to be MCAR if the probability that the responses are missing is independent of both $Y_i^O$ and $Y_i^M$, $P\mathrm{<mo\; stretchy="false">(</mo>}{R}_i\mathrm{<mo>|</mo>}{Y}_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{R}_i\mathrm{<mo\; stretchy="false">)</mo>.}$ The missing data can be missing arbitrary or non-arbitrary. Here for MCAR mechanism, we will set up arbitrary missingness. In the simulation design, we assume the missingness is 25% at time 2 and 25% at time 3. Because we generate drop out missingness for the MAR and MNAR mechanisms, we checked that the drop out patterns in our MCAR data satisfied MCAR. Let $D_i\mathrm{<mo>=</mo>0}$ for subject $i$ who has data at all time points $t\mathrm{<mo>=</mo>1,2,3}$ (no drop out), $D_i\mathrm{<mo>=</mo>1}$ for the dropout case at time 1 and $D_i\mathrm{<mo>=</mo>2}$ when the dropout occurs at time 2. Also, let $last$ denote the last observed value of the response for subject $i$. Then, a logistic regression model for each outcome is done separately as follows:

$log\mathrm{<mo\; stretchy="false">[</mo>}\frac{P\mathrm{<mo\; stretchy="false">(</mo>}{D}_i\mathrm{<mo>=</mo>}j\mathrm{<mo>|</mo>}{D}_i\ge j\mathrm{<mo\; stretchy="false">)</mo>}}{1-P\mathrm{<mo\; stretchy="false">(</mo>}{D}_i\mathrm{<mo>=</mo>}j\mathrm{<mo>|</mo>}{D}_i\ge j\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">]</mo><mo>=</mo>}{\beta }_{0j}+{\beta }_{1}las{t}_i+{\beta }_2{X}_i+{\beta }_{3}las{t}_i{t}_j$

The MCAR mechanism implies that ${\beta }_1$ and ${\beta }_3$are zero. Hence, a test for MCAR is as test of the null hypothesis $H_0\mathrm{<mo>:</mo>}{\beta }_1\mathrm{<mo>=</mo>}{\beta }_3\mathrm{<mo>=</mo>0}$, which is not rejected when ${\widehat{\beta }}_1$ and ${\widehat{\beta }}_3$ are both not statistically different from zero. If the simulated dataset led to rejection of $H_0$, we discarded those data and generated a new data set until $H_0$ was not rejected.

Missing at random MAR

The missing at random mechanism implies the missingness probability is independent of $Y_i^M$, i.e., $P\mathrm{<mo\; stretchy="false">(</mo>}{R}_i\mathrm{<mo>|</mo>}{Y}_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}R\mathrm{<mo>|</mo>}{Y}_i^O\mathrm{<mo\; stretchy="false">)</mo>.}$ Here the design is much easier than for MCAR. We just design the missingness to be related the observed data and independent of unobserved data. In this simulation design, we set up the drop out after the first time point. If the logit of the response is lower than a specified cutpoint, then the subject drops out at the next time point for all the subsequent times.

Missing not at random MNAR

The missing not at random mechanism holds when the probability the responses are missing is not independent of either $Y_i^M$ or $Y_i^O$. Then we set the missingness related to the observed and unobserved data. After the first time point, if the logit of the response is lower than a specified cutpoint, then the subject drops out at that time point for all the subsequent times. Here the cause of the missingness comes from observed and unobserved data, satisfying MNAR. To address the multivariate structure, we further classified mechanisms as “both missingness" when missingness occurs in both outcomes, and “partial missingness" when only one outcome is incomplete. The term “mixed missingness" when the two outcomes have different missing mechanisms. Defining different missing patterns is in order to study the changes in the regression effects estimations and the correlation parameters. In the following table, there are eight missing data mechanisms (Figure 1) (Table 2).

Figure 1 Test power versus sample size.

Handling missing data

There are many traditional methods of handling missing data in longitudinal studies. We conducted four methods as follows.

Completed case analysis

This method is defined by deleting all subjects who have missing data. It includes only the completed cases. The primary condition to use this method for valid inference is the MCAR assumption. Based on Little,¹ It will yield unbiased estimates of the mean responses when the assumption of the missing data is MCAR. However, the completed case analysis is not appealing when the size of the completed subjects is small relative to the whole number of subjects because it will reduce the statistical power and increase the standard error of the estimates Allison.⁸ Also, if the data is not MCAR, then the parameter estimates could be biased.

Imputation

The imputation approach is widely used in practice. The basic idea is filling the missing values with imputed values. The imputed values are chosen by many methods. We will discuss three methods: Mean Substitution (MS), Last Observation Carried Forward (LOCF) and Regression Imputation (RI). The advantage of using imputation methods is the possibility to use the standard statistical analysis methods for the completed data after the imputation.

Mean substitution (MS)

In this method, we easily use single imputation for each missing data by the mean of available observed values. This method may causes bias in the parameter estimations and under estimate the variability Cook et al.¹⁰ Thus, it may artificially reduce the variability of the parameter estimations. In this paper, we accommodate the multivariate structure by imputing the mean of each outcome, for each occasion, and for the same treatment. That means when we have two outcomes, three occasions and treatment binary variable, we will impute the missing values by 12 means instead of single mean. Each missing value is imputed by the mean of its outcome, its occasion, and treatment.

Last observation carried forward (LOCF)

The LOCF is frequently used in clinical trials, it is a single imputation method that is widely used for longitudinal data, despite that it leads to bias expect for MCAR mechanism. In this method, each missing value is replaced by the last observed value for that subject. It is easy to impute LOCF method, but it is unrealistic to assume the observations following the dropouts remain unchanged unless the dropout is due to cure. Shao & Zhong¹¹ conducted studies and proved that LOCF method could have biased results and reduce the precision of the parameter estimations.

Regression imputation (RI)

The idea behind the regression imputation is simple. It is based on assuming the distribution of the missing responses is close to the distribution of the observed responses. Here, we impute the missing responses using the prediction from the regression equation derived from the observed responses. Then, we remodel the completed responses after the imputation, using the prediction method to impute the missing responses seems not appealing when the cause of the drop out is related to the efficiency of the treatment. For example, when the sick subjects tends to drop out more than the healthy subjects, this leads to different estimates of the observed and missing responses. In MAR and MNAR assumptions, the there is a statistical difference between responses means of the observed and missing responses.

Inverse probability weight method (IPW)

The GEE method can be adapted to treat the dropout missingness for MAR mechanism. One applicable method is the inverse probability weighted (IPW) GEE approach of Preisser et al.¹² and Robins & Rotnitzky.¹³ It provides a method to handle the dropout missingness in the lack of good reasons to assume the MCAR mechanism. The idea is to weight the observed data to account for the probability of unobserved data. The weights are obtained by estimating the probability of dropout as a function of observed responses prior to dropouts and covariates. The probability of not dropping out is called the propensity score. Here in the multivariate longitudinal data, the appropriate way is to estimate the propensity score weights for each outcome $k$. Let ${\pi }_{ijk}$ denote the conditional probability of subject $i$ being observed at time $j$ for outcome $k$ given the readings history of that subject at prior times.

${\pi }_{ijk}\mathrm{<mo>=</mo>}P\mathrm{<mo\; stretchy="false">(</mo>}{R}_{ijk}\mathrm{<mo>=</mo>0}\mathrm{<mo>|</mo>}{R}_{i1k}\mathrm{<mo>=</mo>}{R}_{i2k}\mathrm{<mo>=</mo>}\dots \mathrm{<mo>=</mo>}{R}_{i\mathrm{,}j-\mathrm{1,}k}\mathrm{,}{X}_i\mathrm{,}{Y}_{i1k}\mathrm{,}{Y}_{i2k}\mathrm{,}\dots \mathrm{,}{Y}_{i\mathrm{,}j-\mathrm{1,}k}\mathrm{<mo\; stretchy="false">)</mo>}$

The appropriate weight for$Y_{ijk}$ is the inverse of the unconditional probability of $Y_{ijk}$being observed, this probability being the cumulative product of${\pi }_{ijk}$ for $j\mathrm{<mo>=</mo>1,2,..,}J$:

$w_{ijk}\mathrm{<mo>=</mo>}\frac{1}{{\pi }_{i1k}{\pi }_{i2k}\mathrm{...}{\pi }_{iJk}}$

In our simulation context, the model for missingness estimates the probability of not dropping out at given time point, given the previous response, treatment values and their interaction. The logistic regression model for each outcome$k$ is:

$logit\mathrm{<mo\; stretchy="false">(</mo>}{\pi }_{ijk}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\beta }_{0k}+{\beta }_{1k}{t}_{ij}+{\beta }_{3k}{X}_i+{\beta }_{4k}{Y}_{i\mathrm{,}j-\mathrm{1,}k}+{\beta }_{5k}{X}_i{Y}_{i\mathrm{,}j-\mathrm{1,}k}\mathrm{.}$

Thus score propensity is estimated for each subject at each time. Then, we analyze the data using weighted GEE model of Shelton⁷ using the “working independence" correlation matrix.

Simulation evaluation

We used five criteria to measure the performance of each scenario and mechanism for the correlated binary data. These are the average, standard deviation (std), the bias to express the difference between the true and estimated parameters, the root of mean square error (RMSE), and the coverage probability (CP) associated with the usual 95% confidence interval. Using the RMSE is good because it is the square root of MSE that measures the precision of the estimates. The coverage probability is the proportion of the nominal 95% confidence intervals from simulated datasets that contain the true parameter values.

The general strategy of this study is generating the artificial multivariate longitudinal data for two outcomes and three occasions, and five correlation scenarios based on the model in 2.1. Then, we set up the eight missing mechanisms described in 2.3. First, The incomplete datasets are analyzed, and the parameters are estimated. The correlation parameters are estimated using sample correlation values for each mechanism, and for each correlation scenario. In scenario 1, where all correlation parameters are assumed to be zero, the estimated parameters are close to zero.

In scenarios 2 and 3, we induced the correlation between the outcomes in times 1, 2 and 3 respectively,${\alpha }_1$, ${\alpha }_2$and ${\alpha }_3$. Table 3 reports the estimated means of 224 samples of the outcomes parameters. In scenario 2 and scenario 3, the true correlation parameters are 0.6 and 0.9, respectively. The estimated correlation parameters in the completed data exhibit increasing bias over time, and this may be due to existence of the time covariate. Comparing the estimated correlation parameters for the complete and incomplete data for different mechanisms, we notice many points. We found ${\alpha }_1$ did not change for all the missing mechanisms because the baseline is full without missingness. The MNAR assumption has a clear effect on ${\alpha }_2$ and ${\alpha }_3$ over the both, partial and mixed mechanisms. In Table 4, the estimated correlation parameters for time factor ${\upsilon }_1$, ${\upsilon }_2$ and ${\upsilon }_3$ represent the correlation between the responses at (time1, time2), (time1, time3) and (time2, time3) within the first outcome and ${\upsilon }_4$, ${\upsilon }_5$ and ${\upsilon }_6$ within the second outcome. In scenarios 4 and 5, we induced the correlation over the occasions factor. We found the MCAR_MCAR and COMP_MCAR yield results very close to parameters in the completed data but the MCAR_MAR has been affected in scenario 5, the scenario that has strongest correlation over the occasions.

Additionally, the results of the remaining MAR and MNAR for both, partial and mixed missingness are biased and have been affected by the correlation apart from the baseline parameters ${\alpha }_1$ and ${\alpha }_4$. We found the MCAR assumption for both outcomes is the appropriate assumption for the GEE model. Now we can figure out its advantage to keep the correlation parameter close to the correlation parameters in the complete data. Also, we may find the missingness pattern of MAR and MNAR mechanisms could affect the parameter estimations due to the bias in the correlation parameters before we fit the model.

In model 1, there are three regression coefficients associated with the intercept,$X$and time covariates for each outcome after we fit the GEE model assuming the unstructured within subject correlation. The parameter estimations of the X covariate are${\beta }_{11}\mathrm{,}{\beta }_{12}$for scenarios 1, 2 and 4 in Tables 5–7, respectively. Generally from all the scenarios, we found the results of MCAR_MCAR and COMP_MCAR mechanisms are close to the regression coefficients of the completed data without missngness. The MAR_MAR and MNAR_MNAR have shown some bias results. In scenarios 2 and 4, the bias clearly affects the MAR more than MNAR in the parameters of the complete outcome. For the partial mechanisms, we found in scenario 1 that the estimated effects appear unbiased for the for the first outcome, which is complete. We conclude that assuming MCAR, which is implicit for GEE models, is affected by MAR and MNAR when they are mixed, especially with the strong correlation over the occasions. Now, we present the analysis after handling the missing data mechanisms for the incomplete data. We handled the missing data using four different methods: CC, MS, LOCF and RI as they are described in section 2.4. In Table 8 the results of the x covariate parameters estimations after handling the missingness just for the mechanism MCAR_MCAR while the rest of the intercepts and time covariate are in appendix B. Since the appropriate mechanism for the GEE models is the assumption MCAR, we will analyze the handling methods for the mechanisms COMP_MCAR, MCAR_MCAR, COMP_MAR and MCAR_MAR. In Table 8 the results of the estimatedx effects on the log odds of response $P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{ijk}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo>}$for the MCAR_MCAR mechanism. We present four evaluation measurements. We found using the MS and LOCF methods have good means, CP and also less bias in scenario 1 for the estimated parameters of the treatment effects. In scenario 2 where the correlation is induced between the outcomes, we found the methods LOCF and MS have better means close to the true values, good CP and less bias. For the results of scenario 4 where the correlation is induced over the occasions, we found the parameter estimation of the completed data are already biased and we don’t find good method to reduce it over the four methods.

In Table 9 the results of handling the mechanism COMP_MCAR for estimated$X$effects on the log odds of response$P\mathrm{<mo\; stretchy="false">(</mo>}{Y}_{ijk}\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo>}$.This mechanism is a mixture when the first outcome is complete and the second is MCAR. We found the best method in all the scenarios is MS with good means, less standard deviation and bias for the estimated parameters of the two outcomes unless in the parameter${\beta }_{12}$in scenario 4. This estimated parameter is already bias in the completed data and MS method have shown a slight reduction in the bias criteria. We found the estimates of the completed outcome were not affected by the MCAR mechanism in the second outcome regardless of the correlation induced between the outcomes or occasions.

In Table 10 where the mechanism is mixed between the MCAR for the first outcomes and MAR for the second outcome. Here it is a good point to testify the stability of the estimated parameters of the first outcome when the parameters of the second outcome has MAR set up pattern for its missing data. We found the MS method has been good for all the scenarios for the two outcomes with good means and less std. This generalization has exception in the second parameter${\beta }_{12}$in scenario 1 and${\beta }_{11}$in scenario 4 when the there a slight difference between the incomplete data and MS results. Here we indicated the mean substitution imputation reduce the bias of the MAR missing mechanism in the second outcomes and doing good results to impute the first outcome when its mechanism is MCAR. Finally, it seems the MS method has good means and less standard deviation and less bias in most the results of handling the missing data. In Table 11, we have mixed missing mechanisms for the two outcomes. While the first outcome is complete, the second outcome has MAR assumption for its missing data. Here the goal is testify the precision of the parameters of first outcomes if the they will be affected by the MAR assumption, especially when the two outcomes are correlated. Here the the additional method IPW is added since it is appropriate to reduce the bias in the MAR assumption. Generally, we found the RI method has good means, CP and less std for all the scenarios while the CC has the worst results. Based on the bias results we found, RI and MS and IPW reduced the bias from the completed data in scenarios 2 and 4 where correlation is increased. Generally, we can conclude from model parameters the MS and IPW have some good results to treat the missingness in the mechanism COMP_MAR.¹⁴

We presented the analysis of incomplete multivariate binary longitudinal data. The correlation parameters for the complete and incomplete data are estimated. We found the missingness affects the estimated correlation among the occasions more than between the outcomes. Also, we found the estimates of MAR and MNAR of incomplete data are affected by the induced correlation over the occasions and the outcomes. This agrees with results of about the bias of estimates of MAR or MNAR mechanisms for GEE models. After imputing the incomplete data with four missing data handling methods, we conclude many points. The mean substitution based on the multivariate structure mostly has good estimates in the mechanism COMP_MCAR, MCAR_MAR and COMP_MAR. It has been a good remark to find the MS imputation reduced the effects of MAR assumption in the mixed mechanisms and generated mostly less biased results. In the mechanism MCAR_MCAR, the LOCF method has good results. Also, using the weighted GEE for the full and mixed MAR assumptions shows some good results. Generally, we recommend using the mean substitution based on the multivariate structure to impute the missing data. It could be a good future work to use the multiple imputation to handle the missingness in the multivariate structure based on Shelton et al.⁷ model.

None.

Authors declares there is no conflict of interest.

1. Little RJ, Rubin DB. Statistical analysis with missing data. US: John Wiley & Sons; 1987.
2. Myers WR. Handling missing data in clinical trials: an overview. Drug Information Journal. 2000;34(2):525–533.
3. Touloumi G, Babiker A, Pocock S. Impact of missing data due to drop-outs on estimators for rates of change in longitudinal studies: a simulation study. Statistics in medicine. 2001;20(24):3715–3728.
4. Newman DA. Longitudinal modeling with randomly and systematically missing data: A simulation of ad hoc, maximum likelihood, and multiple imputation techniques. Organizational Research Methods. 2003;6(3):328–362.
5. Hening DA. Missing Data Imputation Method Comparison in Ohio University Student Retention Database. USA: Ohio University; 2009.
6. Ali AM, Dawson SJ, Blows FM, et al. Comparison of methods for handling missing data on immuno histochemical markers in survival analysis of breast cancer. Br J Cancer. 2011;104(4):693–699.
7. Shelton BJ, Gilbert GH, Liu B, et al. A SAS macro for the analysis of multivariate longitudinal binary outcomes. Comput Methods Programs Biomed. 2004;76(2):163–175.
8. Allison PD. Multiple regression: A primer. USA: Pine Forge Press; 1999.
9. Alzahrani H. Generating multivariate longitudinal binary random variables for gee models using bridge distribution. International Journal of Advanced Research. 2017;5(2):163–175.
10. Hedeker D, Gibbons RD. Longitudinal data analysis, 451 volume. USA: John Wiley & Sons; 2006.
11. Cook RJ, Zeng L, Yi GY. Marginal analysis of incomplete longitudinal binary data: a cautionary note on locf imputation. Biometrics. 2004;60(3):820–828.
12. Shao J, Zhong B. Last observation carry–forward and last observation analysis. Stat Med. 2003;22(15):2429–2441.
13. Preisser JS, Lohman KK, Rathouz PJ. Performance of weighted estimating equations for longitudinal binary data with drop-outs missing at random. Stat Med. 2002;21(20):3035–3054.
14. Robins JM, Rotnitzky A. Semiparametric efficiency in multivariate regression models with missing data. Journal of the American Statistical Association, 1995;90(429):122–129.