A Copula-based Imputation Model for Missing Data of Mixed Type in Multilevel Data Sets
Jiali Wang Bronwyn Loong Anton H. Westveld Alan H. Welsh


1]Jiali WangCorresponding author: Jiali Wang, Research School of Finance, Actuarial Studies and Statistics, College of Business and Economics Building 26C, The Australian National University, Canberra, ACT 2601, Australia. Phone: +61 2 612 57290. Email: u5298171@anu.edu.au
1]Bronwyn Loong
1,2]Anton H. Westveld
3]Alan H. Welsh


[1]Research School of Finance, Actuarial Studies and Statistics, College of Business and Economics, The Australian national University, AUS 
[2]Statistics Laboratory @ the Bio5 Institute & Statistics G.I.D.P., The University of Arizona, USA
[3]Mathematical Sciences Institute, College of Physical & Mathematical Sciences, The Australian National University, AUS

A Copula-based Imputation Model for Missing Data of Mixed Type in Multilevel Data Sets
    [
    
======================================================================================

[L]December 30, 2023 



We propose a copula based method to handle missing values in multivariate data of mixed types in multilevel data sets. Building upon the extended rank likelihood of <cit.> and the multinomial probit model, our model is a latent variable model which is able to capture the relationship among variables of different types as well as accounting for the clustering structure. We fit the model by approximating the posterior distribution of the parameters and the missing values through a Gibbs sampling scheme. We use the multiple imputation procedure to incorporate the uncertainty due to missing values in the analysis of the data. Our proposed method is evaluated through simulations to compare it with several conventional methods of handling missing data. We also apply our method to a data set from a cluster randomized controlled trial of a multidisciplinary intervention in acute stroke units. We conclude that our proposed copula based imputation model for mixed type variables achieves reasonably good imputation accuracy and recovery of parameters in some models of interest, and that adding random effects enhances performance when the clustering effect is strong.





§ INTRODUCTION

Multivariate analysis often involves understanding the relationship among variables of different types. Our motivating data set is from the Quality in Acute Stroke Care (QASC) study, which implemented a multidisciplinary intervention to manage fever, hyperglycaemia and swallowing dysfunction in acute stroke patients <cit.>. This study was one of the largest rigorously evaluated clinical trials which showed that organised stroke unit care significantly reduced death and disability among stroke patients. There were 19 acute stroke units in New South Wales, Australia that participated in the study, and they were randomly assigned to an intervention group (10 units) and a control group (9 units). A pre-intervention and a post-intervention cohort of patients were recruited , their demographic variables such as age, gender and marital status were obtained, and process of care variables such as temperature, time from onset to hospital and length of stay in hospital were recorded. The researchers were interested to see if the implementation of the protocols reduced death and dependency, and improved physical and mental health scores. The four primary outcome variables considered were: (1) modified Rankin Scale (an ordinal variable ranging from 0 to 6, measuring the degree of disability or dependence in daily activities); (2) Barthel index (an ordinal variable ranging from 0 to 100, which also measures performance in activities of daily living. It is usually reported as a dichotomised variable with 60 or more and 95 or more as cut points); (3) mean SF-36 mental component summary score; (4) mean SF-36 physical component summary score. Mental and physical component summary scores were measured on continuous scales between 0 and 100.

In the QASC study, all the four outcome variables had moderate amounts of missing data and most of the explanatory variables had missing values as well (Table 1). Ignoring all the patients with missing values, which is known as complete case analysis, is a commonly used approach to handle missing data but may lead to biased estimates and reduced statistical power. In other words, the remaining cases may not be representative of the target population if we ignore them completely. The smaller sample size also decreases the power to detect significant treatment effects. Due to the potential for positive dependence among units within the same cluster, this is especially serious in multilevel data sets. Case-wise deletion reduces the sample size of patients within hospitals and the number of hospitals at the same time if any information at the hospital level is missing. As a consequence, both the variations between and within hospitals may not be accurately estimated. An alternative approach  is to `impute' missing values, so that after imputation complete data analysis can be performed using standard software. Some ad-hoc procedures include mean imputation and last observation carried forward.  More principled imputation methods are model-based, such as joint modelling <cit.> and fully conditional specification <cit.>.




Current methods to handle missing data are potentially inadequate to apply to the QASC study which is complicated by the clustering effect and the mix of variable types. <cit.> proposed using a semiparametric copula model based on the extended rank likelihood to analyse multivariate data of mixed types. We extend the work of <cit.> by adding random effects to introduce correlation among individuals within clusters. The model in <cit.> can only be used for continuous and ordinal variables, so we consider a multinomial probit model to handle nominal variables. We then evaluate our model by its ability to recover missing data and estimate the true parameters in some models of interest in both a simulation study and a real data study.

The structure of this manuscript is as follows. In section 2 we briefly summarize some popular multivariate techniques to perform missing data imputation and review the general Gaussian copula model and the extended rank likelihood for semiparametric copula estimation as discussed in <cit.>. In section 3 we describe this extended rank likelihood with random effects and combine the copula model with a multinomial probit model. We outline our algorithm to impute missing data in a multilevel data set using our proposed copula model. In section 4, we present and discuss the results of our simulation and real data studies to evaluate our model. The proposed model is compared against several conventional methods using readily available software packages.  Section 5 provides concluding remarks and discusses some future research.




§ BACKGROUND OF MISSING DATA IMPUTATION

Let Y=(Y_obs,Y_mis) denote the `complete' data, with observed part Y_obs and missing part Y_mis. Let θ denote the parameter describing the `complete' data Y. Throughout this paper we assume the data are Missing at Random (MAR)<cit.>, meaning that the probability of missing an entry only depends on the observed data, not on the entry value itself, so that inference about (Y_mis,θ) can be made based on only the observed data Y_obs, and no extra effort is needed to model the missing data process <cit.>. The MAR assumption cannot be tested except in artificial simulation settings, however, it is a simplifying assumption which can be made more reasonable by expanding the model to include more variables that are related to the missing data. Data augmentation <cit.> is often used as a simulation based computational algorithm to approximate the joint posterior distribution of p(θ,Y_mis|Y_obs). It draws Y_mis from p(Y_mis|Y_obs,θ) and θ from p(θ|Y) iteratively. The θ can be treated as coming from the marginal distribution p(θ|Y_obs) and the Y_mis can be treated as coming from p(Y_mis|Y_obs), if our interest lies in filling in the missing values to create complete data sets.



 §.§ Multiple Imputation

Having obtained guesses for the missing data from an imputation model (which will be discussed further below), we cannot treat them as the `true' data because of the uncertainty due to nonresponse. <cit.> proposed multiple imputation (MI) to obtain M independent draws of Y_mis from p(Y_mis|Y_obs) to create M complete data sets. Combining rules are then applied to the parameter estimates from each of the M complete data sets to obtain a single inferential result, as follows.

Let Q be the target population quantity of interest, for example, the coefficients of a regression model. Suppose q̂_m is the point estimate of Q from the m^th imputed complete data set and ŵ_m is an associated measure of sampling variance, m=1,...,M. Three quantities are required for inference on Q:

    q̅   =1/M∑_m=1^Mq̂_m
    
    B   =1/M-1∑_m=1^M (q̂_m-q̅)^2
    W̅   =1/M∑_m=1^Mŵ_m

The analyst uses q̅ as the point estimate of Q. The sampling variance of q̅ is estimated by T=W̅+ (1+ 1/M) B. The total variance associated with q̅ is a function of the within imputation variance and the between imputation variance.

Next we discuss common approaches to impute missing values.



 §.§ Approaches to Generate Imputations for Missing Values

A good imputation method aims to preserve relationships among survey variables of interest. The joint modelling (JM) approach usually assumes the data follow an elliptical joint distribution, for example, a multivariate normal or a multivariate t distribution. For continuous variables, some transformations may be needed to approximate the assumed distribution <cit.>. Discrete variables are treated as if they were generated from the underlying continuous variables and then discretized. Most software packages implement the joint modelling approach by first transforming any variables with missing values into responses that follow a multivariate normal distribution. The transformed responses are then regressed against the fully observed variables. The software packages that implement this approach include norm <cit.> and Amelia <cit.> in R and PROC MI in SAS. Other joint modelling techniques include loglinear models and general location models specifically designed for categorical data and mixed data respectively <cit.>. Another useful package in R - pan <cit.> is designed to impute missing values in panel data, assuming a multivariate Gaussian distribution with random effects. <cit.> further extended Schafer's multilevel imputation model by allowing for multivariate response variables at all levels of a data hierarchy, and used Box-Cox type normalizing transformations for continuous non-Gaussian responses. Although elliptical distributions allow for parsimonious description of data, they are restrictive in the marginal distributions which are fully determined by the parent joint distribution and are restrictive in capturing complex dependencies among variables.

The fully conditional specification (FCS) <cit.> approach breaks the joint model into a series of univariate regression models. Generalized linear models are often specified to accommodate different types and shapes of variables as well as adding constraints among variables. This method has been implemented by many software packages, for instance, mice <cit.> and mi <cit.> in R, ice in STATA <cit.> and a SAS-based software IVEware <cit.>. To the best of our knowledge, there are no available packages to implement the multilevel fully conditional specification except for the `mice.impute.2l.norm' function in the mice package in R, which fits mixed effects linear regression models for variables with missing values. Because of the lack of packages for practitioners, some authors have investigated including indicator variables for clusters <cit.> in the imputation models or ignore the clustering effects. The main criticism of the fully conditional methods, however, is the lack of theoretical justification to ensure the univariate conditional distributions converge to a proper joint distribution.

Several papers have compared JM and FCS MI, but there is no clear conclusion under which circumstances practitioners should favour one over the other. <cit.> performed simulations under three missing data mechanisms and their results showed that JM and FCS produce similar results despite the data not being multivariate normal. <cit.> not only assessed the accuracy of the coefficients fitted to models of interest, but also the accuracy of imputed values. Their study found that FCS imputed more accurately for categorical variables than JM but the differences were small for continuous variables. <cit.> studied the performance of JM and FCS in multilevel settings, and showed using simulations that FCS MI outperforms JM MI in having less bias, and when the intraclass corrlation is small, more accurate parameter estimates are obtained from both JM and FCS.



 §.§ Copulas

To provide more flexibility in the marginal distributions while at the same time ensuring a proper joint distribution, we consider copula modelling approaches to impute missing values. The word `copula' means `a link, tie, bond'. In mathematics and statistics, it means joining together one-dimensional distribution functions to form a multivariate distribution function. Specifically, the distribution functions for the random variables y_1,...,y_p are F_1(y_1),...,F_p(y_p). Sklar's theorem <cit.> shows that there always exists a function C, such that, F(y_1,...,y_p)=C(F_1(y_1),...,F_p(y_p)), where the function C is called the copula function. Each of the variables is modeled by the marginal distribution F_l(y_l)=u_l, l=1,...,p, which is uniformly distributed, and their dependence is captured by the copula function C. Copula modelling has proven to be very powerful for modeling variables of different types and shapes, when there is an underlying dependence among them. It adopts a `bottom-up' strategy where the starting point is the marginal distributions F_l, which are then glued together by the copula function C. In the `top down' joint modelling approach, the marginal distributions are fully determined by their parental joint distribution so that there is no flexibility to model them. In addition, copula models guarantee the existence of a compatible joint distribution which is not guaranteed by the fully conditional specification approach. Existing models, like multinomial (ordered) probit models for (ordered) categorical data can be treated as special cases of copulas, because the underlying latent variables corresponding to each category are assumed to follow a multivariate Gaussian distribution <cit.>.

In a copula model, the parameters are the marginal distributions F_l and the copula function C. <cit.> developed a fully Bayesian estimation procedure to model the joint distribution of both sources of parameters. However, specifying each of the marginal distributions is labour intensive and variables in real data sets may not be accurately represented without a large number of parameters. Some authors suggested transforming the variables using the empirical distribution F̂_l to get pseudo data <cit.> and avoid the parametric estimation of marginal distributions. However, this only applies to continuous variables. To link the discrete variables with continuous latent variables, <cit.> provided a simple way of analysing the correlation among variables with meaningful ordering (continuous and ordered categorical variables), via the extended rank likelihood. This makes use of the fact that the order of the underlying latent variable is consistent with the observed data, and inference about the association parameters can be drawn from the `rank-based' latent variables through a simple parametric form. The extended rank likelihood has already been applied to other closely related models, for example, a general Bayesian Gaussian copula factor model proposed by <cit.> and a bifactor model considered by <cit.>, can be treated as imposing a special structure on the correlation matrix of a Gaussian copula.

Using the copula model as an imputation engine is relatively new but has drawn some attention in the literature. <cit.> were among the first authors to consider imputation using a Gaussian copula where the missing data pattern was monotone. <cit.> found that copula based imputation from the Archimedian family compared favourably with nearest neighbour donor imputation and regression imputation by the EM algorithm. <cit.> compared the performance of imputation by the copula model using the extended rank likelihood approach <cit.> with JM (as implemented in Amelia) and FCS (as implemented in MICE) and concluded that the copula imputation approach maintained the prediction accuracy at least as well as the other two approaches but with faster convergence of the sampling algorithm.




§ SEMI-PARAMETRIC GAUSSIAN COPULA MODEL



 §.§ The Extended Rank Likelihood with Random Effects

Among a variety of copulas, we focus on the Gaussian copula in this paper. For further theoretical details of copulas, see <cit.> and for a good summary of some applications of copulas, see <cit.>.  Rather than assuming a Gaussian distribution on the data Y directly, the Gaussian copula specifies a joint multivariate Gaussian distribution on the corresponding latent variables as defined next. Let l=1,...,p denote the index of the l^th random variable. Then the l^th latent variable is z_l=Φ^-1(u_l), where u_l=F_l (y_l). That is, C(u_1,...,u_p|Γ)=Φ_p(Φ^-1(u_1),...,Φ^-1(u_p)|Γ)=Φ_p(z_1,...,z_p|Γ), where Φ_p(·|Γ) is the cumulative distribution function of the p-variate normal distribution, with mean zero and correlation matrix Γ. Note that the Gaussian copula can reach the full range of pairwise correlation (-1,1) and the parameters that need to be estimated only come from the correlation matrix Γ.

<cit.> derived a rank-based likelihood to estimate the correlation matrix Γ so that there is no need to specify the marginal distributions F_l. The idea is that since we know Φ^-1(F(·)) is a monotone transformation, the ordering of data Y provides partial information about what z should be, that is, y_i_1l<y_i_2l implies z_i_1l<z_i_2l. Suppose we have in total N observations, n=1,...,N. Observing y=(y_1,...,y_N) tells us that z=(z_1,...,z_N) must lie in the set: {z ∈ℝ ^N × p: max{z_hl: y_hl<y_nl}< z_nl < min {z_hl:y_hl> y_nl}}. Let `D' denote the set of all possible z which is consistent with the ordering of y. Then the event `z ∈ D' can be treated as the observed event upon which inference of Γ is made. The full likelihood can be decomposed as

    p(y|Γ, F_1,...,F_p)   =p(z∈ D, y|Γ, F_1,...,F_p)
       =p(z∈ D|Γ) × p(y|z∈ D,Γ,F_1,...,F_p).


<cit.> proved that it is partial sufficient (in the sense of G-sufficient and L-sufficient) to carry out inference about Γ based on the density p(Z∈ D|Γ) and he referred to it as the `extended rank likelihood'. In doing so, we lose the information about Γ from the density p(y|z∈ D,Γ,F_1,...,F_p), but we do not need to estimate the potentially complicated marginal distribution functions and the extended rank likelihood provides a more general and flexible framework for joint modelling.


To take into account clustering effects, we extend Hoff's work by adding random effects to the Gaussian copula model at the latent variable level. The idea is that the clustering of the observed data is carried through to the latent variable level.  Our model can be described as

    z_ij|b_i1∼ N_p(b_i1,Γ_1), b_i1∼ N_p(0,Ψ_1),

where i={1,...,m} is the group index, j={1,...,n_i} is the individual index within group i, Γ_1 is a correlation matrix and Ψ_1 is a variance-covariance matrix for z_ij and b_i1. Both z_ij and b_i1 are vectors of length p, because we are considering l={1,...,p} variables jointly. In this model, the parameters that need to be estimated are in (Γ_1,Ψ_1), which can be thought of as splitting the total correlation into two parts, the variability within groups and the variability between groups. However, like any model that relies on the ordering of the data but not their magnitude,  model (<ref>) suffers from an identifiability problem without constraints on Γ_1. To see this, if we shift the location of the latent variable z_l by μ_l and scale it by σ_l, the model remains unchanged because the new latent variables satisfy the order of the observed data as well. The extended rank likelihood contains only the information about the relative ordering of z but no information about their location and scale. To solve the identifiability problem of scale, we fix Γ_1 to be a correlation matrix instead of a covariance matrix. In other words, there is no need to estimate the variances of z conditional on the random effects, so we fix them as 1. Because the marginal distributions of z have mean equal to 0, there is no identifiability issue for location. We will briefly describe how to add covariates in the discussion section so that the mean of z is no longer 0.



 §.§ Copula Model for Mixed Type Variables

Notice that the extended rank likelihood described above only applies to continuous and ordinal variables, since it makes no sense to consider meaningful numeric values for nominal variables (categorical variables without ordering). To include nominal variables in the copula model as well, we consider a multinomial probit model <cit.> which can be treated as a Gaussian copula. The idea is to relate a nominal variable to a vector of latent variables which can be thought of as the unnormalized probabilities of choosing each of the categories. Suppose a single nominal variable y has K categories, and we define K-1 latent variables for unit i as w_i=(w_i1,...,w_i,K-1) which follow a multivariate Gaussian distribution. Since all the variables appear on one side and we model them jointly, there are no covariates as predictors for now, therefore we only include the intercept term β vector to represent the relative differences between each category 1,...,K-1 compared with the baseline category K. To add a second level to the hierarchy, again we have the random effects b_i2 in the model

    w_ij=   β+b_i2+ϵ_ij
    
      b_i2∼    N_K-1(0,Ψ_2),  ϵ_ij∼ N_K-1(0,Γ_2)
    
      y_ij =   k  if  w_ijk >w_ijk' and  w_ijk>0, for k'≠ k 
    
          K  if  w_ijk<0,  for all k=1,...,K-1.

The rule of deciding the category is a mapping from the latent variables vector to the observed category. The category k=1,...,K-1 is observed if the k^th element of the vector w_i is the largest and greater than 0; the last category K is observed if the largest element in w_i is smaller than 0. We also fix the diagonal elements of Γ_2 equal 1 to be identifiable.

To provide a unified framework of multivariate analysis for mixed type variables, we combine model (<ref>) for variables with ordering and model (<ref>) for variables without ordering as follows


    z_ij|b_i1∼ N_p(b_i1,Γ_1), w_ij∼ N_K-1(β+b_i2,Γ_2),
       b_i=(b_i1,b_i2)∼ N_p+K-1(0,Ψ), Ψ=
       [  Ψ_1 Ψ_12; Ψ_21  Ψ_2 ],
       (z_ij,w_ij)|b_i∼ N_p+K-1((0,β)+b_i,Γ), Γ=
      [  Γ_1 Γ_12; Γ_21  Γ_2 ].

The correlations between variables y_1,...,y_p and y_p+1 are modelled through the off-diagonal matrices Ψ_12 and Γ_12 at the group level and the individual level respectively. Since both Γ_1 and Γ_2 have diagonal elements fixed to be 1, the big matrix Γ is an identifiable correlation matrix.



 §.§ A Gibbs Sampler

 A Gibbs sampling scheme is constructed to approximate the joint posterior distribution p(β,Ψ,Γ,b,z,w,y_mis|y_obs) where the unknown quantities in model (<ref>) are the parameters (β,Ψ,Γ) and the latent variables (b,z,w) as well as missing data y_mis. A simple conjugate prior does not exist for a correlation matrix, and we follow the idea in <cit.> of employing a parameter expansion approach <cit.> to facilitate the MCMC sampling. Specifically, we put an Inverse Wishart prior on the matrix Γ̃ which is the semi-conjugate prior in a multivariate Gaussian sampling model. Then the full conditional distribution of Γ̃ can be derived analytically. After updating Γ̃ in each iteration, we rescale it to be a correlation matrix Γ. For ease of computation, we put an improper flat prior on β and a semi-conjugate Inverse Wishart prior on Ψ, where the hyperparameters are the degrees of freedom ν and the scale matrix Λ

    p(β)∝ 1, 
       Ψ∼ Inv Wishart (ν_1,Λ_1),
       Γ̃∼ Inv Wishart (ν_2,Λ_2).


Under these priors, it is straightforward to derive the full conditional distributions for the parameters (β,Γ,Ψ) as follows

  
  * p(β|…) ∼ N(1/N∑_i=1^m ∑_j=1^n_i(w_ij-b_2i-Γ_21Γ_1^-1(z_ij-b_i1)),1/N(Γ_2-Γ_21Γ_1^-1Γ_12));
  
  * p(Γ̃|…) ∼ Inv Wishart (ν_1-1+N,Λ_1+∑_i=1^mϵ_i^Tϵ_i), where ϵ_i=(z_i,w_i)-(0,β)-b_i

  Γ_[g,h]=Γ̃_[g,h]/√(Γ̃_[g,g]Γ̃_[h,h]), g,h=1,...,p, Γ is rescaled from Γ̃ after each sampling;
  
  * p(Ψ|…)∼ Inv Wishart (ν_2+m,Λ_2+B^T B), where B=(b_1,...,b_m).

From the joint Gaussian distribution of (z,w), we can derive the following conditional distributions for the latent variables z,w,b:

  
  * p(z_ij|…) ∼ N(b_1i+Γ_12Γ_2^-1(w_ij-β-b_i2),Γ_1-Γ_12Γ_2^-1Γ_21);
  
  * p(w_ij|…) ∼ N(β+b_2i+Γ_21Γ_1^-1(z_ij-b_i1),Γ_2-Γ_21Γ_1^-1Γ_12);

  z_ij and w_ij should be sampled from a truncated Gaussian distribution and a Gaussian distribution under the observed category constraint respectively, see below for details.
  
  * p(b_i|…) ∼ N (U_i(Γ^-1⊗ 1_n_i^T)vec((z_i,w_i)-(0,β)),U_i),where U_i=(Ψ^-1+n_iΓ^-1)^-1.


The operator ⊗ is the Kronecker product and vec() is the operator that vectorizes a matrix by stacking its columns. Updating the latent variable z is achieved by sampling from a truncated multivariate Gaussian distribution, where the lower and upper bounds for each single entry z_ijl are determined by: lw=max(z_hl:y_hl<y_ijl) and up=min(z_hj:y_hj>y_ijl) respectively, and h is the index that searches over all the rows in the l^th variable. For example, the lower bound for z_ijl is the maximum value of the latent variable z in the l^th column whose corresponding y is smaller than y_ijl and the upper bound can be defined accordingly. Updating the latent variable w is achieved by sampling from a multivariate Gaussian distribution under the constraint of the observed category by an acceptance and rejection algorithm <cit.>. Specifically, we sample a w vector from the multivariate Gaussian distribution and accept this draw if and only if the maximum element of w occurs at the place of the observed category and is greater than 0, or all the elements in w are smaller than 0 and we observe the reference category K. We continue to sample w until a draw is accepted. When there are missing values in (y_1,...,y_p), the lower and upper bounds in z are undefined, and/or any missing value occurs in y_p+1, the observed category in y_p+1 does not exist. In these cases, we just sample z and/or w from the multivariate Gaussian distributions without the constraints.

To sample missing values for variables with ordering, we apply the monotone transformation on z: y_ijl=F̂_l^-1[Φ(z_ijl)], l=1,...,p, where F̂_l is the univariate empirical distribution function of variable y_l. To sample the missing values in nominal variables, we choose the category corresponding to the largest element in w if it is greater than 0, and choose the reference category if the largest element in w is smaller than 0.



§ SIMULATIONS AND REAL DATA ANALYSIS ON THE QASC

We evaluated the performance of the proposed model through two simulation studies: (i) simulated artificial data with missing values and (ii) the QASC data set with randomly deleted records. We compared the proposed imputation model with other commonly used procedures to treat missing data.



 §.§ Simulation Based on Artificial Data

We generated 100 complete multilevel data sets with correlated variables of different types, and then deleted some entries under the MAR assumption. The total number of clusters in each data set was 20, the cluster size was 50, and the five variables X_1, X_2, X_3, X_4, X_5 had Gamma, binary, nominal, ordinal, and normal distributions respectively. The variable X_1 followed a skewed Gamma distribution: X_1 ∼ Gamma(3,0.5). We assumed all the subsequent variables were generated depending on the previous ones, to introduce correlation among variables. Specifically, X_2 was a binary variable such that logit(p_X_2)=X_1+ϵ_ij, where p_X_2 is the probability that X_2 equals 1 and ϵ_ij∼ N(0,1). The nominal variable X_3 had 4 categories and was generated by a multinomial probit model, so that 3 latent variables were needed: (l_X_3,1,l_X_3,2,l_X_3,3) ∼ N((X_1,x_2)B_X_3, C_X_3), where B_X_3 is a randomly generated coefficient matrix of dimension 2×3 and C_X_3 is a correlation matrix of dimension 3 × 3. The category in X_3 was chosen to be k (for k=1,2,3) if l_X_3,k was the largest component and was greater than 0; and was chosen to be 4 if max(l_X_3)<0. Because we aimed to create a data set with a multilevel structure, we let the ordinal variable X_4 be generated from a random intercept model, l_X_4=b_X_4,i+X_1+X_2+ β_X_3 X_3+ϵ_ij, with ϵ_ij∼ N(0,1), b_X_4,i∼ N(0,ρ_X_4), and β_X_3 a vector of length 3, corresponding to the 3 categories in X_3. Three thresholds were used to determine four levels, they were the 20%, 30%, 50% quantiles of l_X_4. Lastly, the normally distributed variable X_5 was also generated from a random intercept model, X_5=b_X_5,i+X_1+X_2+ β_X_3 X_3+β_X_4 x_4+ϵ_ij, where ϵ_ij∼ N(0,1), b_X_5,i∼ N(0,ρ_X_5), and β_X_3 and β_X_4 are vectors of length 3.

To create missing data under the MAR assumption, we assumed X_5 was completely observed and that the probabilities of missingness in X_j (j=1,...,4) depended on X_5. Specifically, let p_mis,ij be the probability that observation i is missing its value for the X_j variable and we assumed that logit(p_mis,ij)=α_jX_5. By adjusting the parameters α_j, we can control the missingness in each variable.

We varied the parameters that generated the data to consider different scenarios: (1) missing rates for each variable from low (10%), median (30%) to high (50%); (2) clustering effect from low (ρ_X_4=ρ_X_5=0.2) to high (ρ_X_4=ρ_X_5=1), corresponding to intra-class correlation coefficients of 0.17 and 0.5 respectively. In the imputation step, we set the number of imputations to be M=10 <cit.>.





 §.§ Simulation Results Summary

To compare the performance of the proposed method with others, we considered six competing methods, some of which have already been implemented in some software packages. These methods are listed in Table 2. We used the package mitools in R <cit.> to implement the combining rules (<ref>) after M complete data sets had been generated.



The assessment of the relative performance of each method was based on the comparison of the imputation accuracy as well as the 95% coverage rates of the coefficients in the following random intercept logistic regression as a model of interest. We chose this model arbitrarily, and believe that the results would also hold for other models of interest.


For each of the 100 simulated complete data sets, we fitted the model logit(p(X_2=1))=b_i+ β_0+ β_1 X_1+ β_2 X_3,2+ β_3 X_3,3+ β_4 X_3,4, b_i ∼ N(0,σ^2). We used the glmer() function in the lme4 package in R to obtain parameter estimates for β=(β_0,β_1,β_2,β_3,β_4). These are our `true' parameter values. After deletion of records by MAR, we applied each of the missing data methods listed in Table 2, and calculated point and variance estimates for β, using the combining rules. We reported the average of the squared bias of the coefficient estimates over the 100 data sets as well as the coverage rates of 95% confidence intervals.

901
0.72

2|c|2*ICC=0.17     3l|CC                                                               3l|JM                                                           3l|FCS                                                               3l|Cluster JM                                                              3l|Cluster FCS                                                          3l|Copula_Hoff                                                           3l|Cluster Copula                                                      
 3-23
2|c|                              10%     30%     1l|50%     10%     30%     1l|50%    10%     30%     1l|50%     10%     30%     1l|50%     10%     30%     1l|50%     10%     30%     1l|50%     10%     30%     1l|50% 
 5*Bias            1l|β_0           0.033                         0.146                         1l|0.552                         0.037                         0.129                         1l|0.388                         0.014                         0.067                         1l|0.198                         0.041                         0.125                         1l|0.339                         0.011                         0.040                         1l|0.096                         0.014                         0.055                         1l|0.104                        0.013                         0.054                         1l|0.092                     

                                 1l|β_1            0.040                         0.183                         1l|0.309                         0.044                         0.133                         1l|0.308                         0.024                         0.092                         1l|0.186                         0.022                         0.088                         1l|0.399                         0.026                         0.105                         1l|0.244                         0.029                         0.104                         1l|0.316                         0.026                         0.095                         1l|0.187                     

                                 1l|β_2            0.034                         0.076                         1l|0.190                         0.026                         0.050                         1l|0.104                         0.022                         0.055                         1l|0.106                         0.033                         0.055                         1l|0.183                         0.021                         0.043                         1l|0.099                         0.024                         0.053                         1l|0.080                         0.021                         0.049                         1l|0.086                     

                                 1l|β_3            0.018                         0.073                         1l|0.164                         0.016                         0.048                         1l|0.155                         0.010                         0.056                         1l|0.118                         0.021                         0.051                         1l|0.140                         0.011                         0.038                         1l|0.079                         0.012                         0.046                         1l|0.080                         0.013                         0.051                         1l|0.087                     

                                 1l|β_4            0.045                         0.126                         1l|0.373                         0.032                         0.194                         1l|0.442                         0.017                         0.098                         1l|0.261                         0.020                         0.089                         1l|0.383                         0.011                         0.047                         1l|0.106                         0.013                         0.057                         1l|0.131                         0.012                         0.055                         1l|0.122                     
 1-1
1-1
5*Coverage     1l|β_0     90      77     1l|67     87     78     1l|79     90      83      1l|83     89      76     1l|73     93      89     1l|87     90      86      1l|80     92      90      1l|87 

                              1l|β_1     97      84     1l|77     95     90     1l|79     98     92      1l|81     98     94     1l|77     94     90     1l|81     89      85      1l|79     100     96      1l|84 

                              1l|β_2     100     92     1l|87     98     92     1l|91     95      100     1l|96     97      89     1l|88     100     98     1l|98     100     100     1l|99     99      100     1l|98 

                              1l|β_3     95      94     1l|82     93     94     1l|88     93      95      1l|88     91      90     1l|76     91      91     1l|85     94      93      1l|87     91      92      1l|92 

                              1l|β_4     87      80     1l|63     89     76     1l|69     93      90      1l|83     90      82     1l|79     90      89     1l|87     95      87      1l|89     93      88      1l|85 
 1-1

tableA comparison of squared bias and coverage of coefficient estimates of the model of interest under seven methods to handle missing data, with ICC=0.17 and missing rates=10%, 30% and 50%.



901
0.72

2|c|2*ICC=0.5     3l|CC                                                               3l|JM                                                           3l|FCS                                                               3l|Cluster JM                                                              3l|Cluster FCS                                                          3l|Copula Hoff                                                           3l|Cluster Copula                                                      
 3-23
2|c|                             10%       30%       1l|50%      10%        30%        1l|50%        10%       30%       1l|50%      10%       30%       1l|50%       10%         30%        1l|50%        10%        30%        1l|50%        10%         30%        1l|50%    
 5*Bias            1l|β_0            0.051      0.164      1l|0.539     0.056       0.160       1l|0.481       0.029      0.114      1l|0.216     0.053      0.197      1l|0.411      0.025        0.074       1l|0.129       0.030       0.100       1l|0.139       0.024        0.076       1l|0.127   

                                 1l|β_1            0.035      0.210      1l|0.542     0.049       0.143       1l|0.424       0.019      0.084      1l|0.212     0.029      0.102      1l|0.373      0.022        0.136       1l|0.303       0.036       0.191       1l|0.364       0.031        0.183       1l|0.290   

                                 1l|β_2            0.040      0.098      1l|0.235     0.068       0.180       1l|0.327       0.027      0.070      1l|0.122     0.034      0.173      1l|0.391      0.023        0.069       1l|0.111       0.028       0.077       1l|0.093       0.022        0.060       1l|0.094   

                                 1l|β_3            0.035      0.105      1l|0.260     0.115       0.250       1l|0.250       0.024      0.088      1l|0.162     0.035      0.119      1l|0.221      0.024        0.071       1l|0.108       0.023       0.073       1l|0.127       0.024        0.071       1l|0.103   

                                 1l|β_4            0.089      0.303      1l|0.492     0.098       0.319       1l|0.514       0.044      0.200      1l|0.340     0.045      0.196      1l|0.456      0.038        0.138       1l|0.175       0.048       0.178       1l|0.211       0.040        0.138       1l|0.183   

1-1
5*Coverage     1l|β_0     91      76     1l|74     87     80     1l|78     89      80     1l|79     88      89     1l|90     90      85     1l|88     90      89     1l|85     92      91     1l|91 

                              1l|β_1     92      80     1l|65     94     86     1l|79     93      92     1l|78     96      92     1l|73     95      90     1l|74     88      83     1l|78     95      93     1l|80 

                              1l|β_2    100     92     1l|80     97     91     1l|80     95      96     1l|92     94      90     1l|86     98     98     1l|93     100     97     1l|97     100     99     1l|91 

                              1l|β_3     89      89     1l|82     90     76     1l|74     90     95     1l|86     90     88     1l|78     88      87     1l|87     90      87     1l|87     89      92     1l|90 

                              1l|β_4     88      76     1l|66     87     81     1l|70     92      84     1l|73     90     83     1l|78     91      88     1l|82     90      82     1l|84     90      84     1l|82 
 1-1

tableA comparison of squared bias and coverage of coefficient estimates of the model of interest under seven methods to handle missing data, with ICC=0.5 and missing rates=10%, 30% and 50%.



Table 3 summarizes the results of the simulation experiments under the three missingness rates (10%, 30% and 50%) using the seven methods, when the ICC used to generate the variables X_4 and X_5 is 0.17. When the missingness rate is 10%, all the approaches give reasonably good results in terms of achieving the nominal coverage rate - 95%, though CC and the two joint modelling approaches (JM and Cluster JM) do worse than the others. A possible reason for this is that in the joint modelling approaches, multivariate Gaussian distributions were specified and this is clearly not true in our data generating process, whereas in the sequential imputation approaches (FCS and Cluster FCS) more flexible univariate imputation models were allowed to best accommodate different variable types. For the copula-based methods, the empirical distribution function transformations were applied before fitting a multivariate Gaussian distribution on the latent variable scale where the dependence among the variables was captured. In addition, the squared bias increases with an increase in missingness rate as expected. With a moderate to high level of missingness, Cluster FCS and our proposed method (Cluster Copula) tend to outperform FCS and Copula_Hoff. While all the methods suffer from under-coverage when the missing rates are 30% and 50%, CC seems to be the worst, producing the most biased results. The results meet our expectation because as the percentage of missing data increases, there is less observed data available to capture the complex dependency among variables. Under the MAR assumption, CC causes the most biased results by only using the complete records while its alternatives make use of all the observed data.

Table 4 is similar to Table 3 except that the performance is evaluated at ICC=0.5. In other words, the data sets exhibit higher levels of clustering. Compared to the results in Table 3, the results are worse across all methods for the higher ICC value. The imputation methods which take into account clustering effects almost always do better than their counterparts, which is not that obvious in Table 3 when ICC=0.17. Conditional imputation methods do better than joint modelling approaches, and the two copula-based methods tend to achieve the best results, for almost half of the simulation settings with the smallest squared bias.

We also compared the imputation accuracy. That is, for each data value we calculated the discrepancy between the average of the 10 imputed values and the before-deletion true values. Note that this comparison is not applied to the CC method.  The Euclidean distance was used to measure the imputation accuracy in the continuous variable X_1 and the ordinal variable X_3: 1/10∑_m=1^10∑_i=1^N (X_i,true-X_i,imp^(m))^2/#miss X, and the misclassification rate was used to measure the imputation accuracy in the binary variable X_2 and the nominal variable X_4: 1-1/10∑_m=1^10∑_i=1^N 1_(X_i,true=X_i,imp^(m))/#miss X.

Figure 1 shows the results of the imputation accuracy for each simulation study. The points are the means of the Euclidian distances/misclassification rates over the missing observations in a single data set, and the error bars show the 5% and 95% quantiles over the 100 data sets. For variable X_1 which follows a Gamma distribution, there is not much difference in imputation accuracy over the six methods. For the nominal variable X_3 our proposed Cluster Copula method always performs the best except for the top-left panel, while the JM approach is always the worst. The misclassification rates for the binary variable X_2 are smallest in all the scenarios when using our proposed Cluster Copula model but do not differ much from those of the other methods. The misclassification rates for the ordinal variable X_4 are again highest for JM and the rates for the copula-based methods are smaller than the others when the missing rates are 30% and 50%. Generally speaking, the copula based methods tend to impute more accurately for categorical variables but also do no worse than other methods for continuous variables. The joint modelling methods, especially JM, give the least accurate imputation as the multivariate Gaussian distribution assumption does not hold. As the missingness rate and/or ICC increase, all the methods for every variable perform comparatively worse in terms of having a larger disparity compared with the true values and higher misclassification rates, but the patterns of relative performance between the six methods remains broadly the same.





 §.§ Simulation Based on the QASC Data Set

We also ran simulation studies using the QASC data set to evaluate our proposed method and other competing methods. Here we treated all the complete cases in the QASC data set (75.34% of the original data set) as the `true' data, and sub-sampled 300 patients, 100 times to create 100 sub data sets. Then for each of the sub data sets, missing values were created, trying to mimic the missing data pattern in the original data set. We distinguish between the demographic variables which we treat as MCAR and the process of care variables which we treat as MAR. Specifically, for the demographic variables: `ATSI', `age', `education' and `marital status', values were randomly deleted to roughly match the missingness percentages in Table 1. For the process of care variables and outcome variables, we assumed their missingness depended on the completely observed variables. A missing indicator variable was associated with every variable with missing data which equaled 1 if an entry was missing. For the missing indicators, we fitted logistic regression models on the original data set for `time taken to hospital', `mean temperature', `modified Rankin Scale', `Bartell Index', `physical health score' and `mental health score' respectively against `gender', `period' and `treatment', and the probabilities of missingness for the sub-sampled data sets were decided by the predicted values of these logistic regression models. We noticed that 9.39% of `Bartell Index', `physical health score' and `mental health score' were missing together, and we also took this into account when creating missing data.

The relative performance for each method was also compared based on the average imputation accuracy and the squared bias and 95% coverage rate of interval estimates of parameters for some models of interest. Ten imputations were created for all the six imputation methods. The accuracy is shown in Table 5. All the discrepancies between the imputed values and the true values were measured by Euclidean distance except for the nominal variables `marital status' and `ATSI' which used the misclassification rates. Our proposed imputation model achieves the smallest disparity more than half of the time (7/11) and Copula_Hoff is superior in performance to the other four methods. It is interesting to note that joint modelling methods perform better than their FCS counterparts (JM vs.FCS and Cluster JM vs.Cluster FCS) and adding clustering effects enhances the imputation accuracy.



The models of interest are based on the models fitted in <cit.>. They fitted logistic regression models for the dichotomous outcomes - `Bartell Index' with cut points equal to 60 and 95, and `modified Rankin Scale' with cut point equaled to 2; and linear models for the continuous variables `physical health score' and `mental health score', including as predictors the variables `treatment', `period' and the interaction between `treatment' and `period'. The models are


    (mrs2)=b_i+β_0+β_1 period+ β_2 treatment +β_3 treatment*period,
       (bi60)=b_i+β_0+β_1 period+ β_2 treatment +β_3 treatment*period,
       (bi90)=b_i+β_0+β_1 period+ β_2 treatment +β_3 treatment*period,
       mcs=b_i+β_0+β_1 period+ β_2 treatment +β_3 treatment*period+ϵ,
       pcs=b_i+β_0+β_1 period+ β_2 treatment +β_3 treatment*period+ϵ.

The coefficient β_3 and its p-value were used to see if the pre-post change in the intervention group was statistically significant to the change in the control group.  All the models included a random intercept term, b_i, to capture the clustering effects.


We first fitted the five models of interest on the completely observed patients in each of the 100 sub data sets, and obtained the parameter estimates β=(β_0,β_1,β_2,β_3) and treated them as the true values. Then the parameter estimates from all the seven competing methods were compared against the true parameters, and the 95% coverage rates were obtained from the 100 repetitions. The results are reported in Table 6. The CC approach has the largest bias and the smallest coverage rate. This is not unexpected because the missing data were generated under the MAR assumption and by CC only about 40% of the data were used to fit the models so that the coefficient estimates are biased with large uncertainty. The proposed method Cluster Copula and Copula_Hoff outperform the other methods with Copula_Hoff doing marginally better than Cluster Copula for the first and second logistic models `mrs2' and `bi60', and Cluster Copula doing better for the fifth linear model for `pcs'. There is little difference between the two copula based methods, because the clustering effects were small in the QASC data set (ICC in the models of interest lay between 0.009 and 0.026), and only one nominal variable (marital status) was considered in the imputation models but did not enter into the models of interest later. In other words, taking the clustering effect into account and giving special treatment to the nominal variable does not affect the inference too much in this case. However, we do observe that when ICC is higher in the variable `pcs', our proposed model achieves better imputation accuracy.

901
		0.6
			
				2|l|2*     3l|CC    3l|JM    3l|FCS    3l|Cluster JM     3l|Cluster FCS                                                            3l|Copula_Hoff    3l|Cluster Copula  
 3-23
				2|l|                      Bias     SD     1l|Coverage     Bias     sd     1l|Coverage     Bias     SD     1l|Coverage     Bias     SD     1l|Coverage     Bias     sd     1l|Coverage     Bias     SD     1l|Coverage     Bias     SD     1l|Coverage 
 4*Modified Rankin Scale 2        1l|β_0       0.029     0.426     1l|91     0.007     0.282     1l|91      0.013     0.284     1l|92     0.006     0.283     1l|90     0.01      0.274     1l|92      0.006     0.282     1l|92      0.007     0.285     1l|92 

				1|l|     1l|β_1   0.049     0.671     1l|94     0.017     0.43      1l|96      0.021     0.431     1l|96     0.015     0.433     1l|98     0.019     0.409     1l|94      0.013     0.432     1l|98      0.016     0.434     1l|97 

				1|l|     1l|β_2   0.022     0.534     1l|93     0.013     0.344     1l|100     0.014     0.344     1l|99     0.014     0.343     1l|98     0.014     0.391     1l|100     0.011     0.345     1l|100     0.013     0.344     1l|99 

				1|l|     1l|β_3   0.043     0.827     1l|88     0.026     0.521     1l|91      0.033     0.523     1l|89     0.023     0.52      1l|90     0.032     0.524     1l|90      0.022     0.522     1l|93      0.024     0.518     1l|92
1-24*Bartell Index 60        1l|β_0        0.051     0.584     1l|87     0.043     0.315     1l|89     0.032     0.329     1l|90     0.048     0.316     1l|86     0.035     0.315     1l|88     0.012     0.327     1l|91     0.021     0.333     1l|90 

				1|l|     1l|β_1   0.028     0.586     1l|97     0.027     0.478     1l|95     0.024     0.487     1l|97     0.022     0.480     1l|97     0.024     0.495     1l|95     0.020     0.495     1l|97     0.026     0.501     1l|95 

				1|l|     1l|β_2   0.043     0.819     1l|85     0.035     0.395     1l|89     0.029     0.405     1l|92     0.025     0.396     1l|93     0.030     0.398     1l|90     0.016     0.411     1l|93     0.018     0.415     1l|93 

				1|l|     1l|β_3   0.083     0.746     1l|90     0.055     0.606     1l|92     0.072     0.624     1l|87     0.046     0.610     1l|95     0.048     0.635     1l|97     0.038     0.634     1l|98     0.055     0.643    1l| 98 
1-24*Bartell Index 95        1l|β_0        1.228     1.293     1l|74     0.148     0.881     1l|81     0.163     0.874     1l|80     0.145     0.875     1l|80     0.160     0.792     1l|78     0.117     0.608     1l|85     0.123     0.627     1l|82 

				1|l|     1l|β_1   3.252     1.379     1l|90     2.518     0.752     1l|96     3.061     0.654     1l|98     1.945     0.733     1l|96     1.854     0.780     1l|95     1.770     0.627     1l|97     1.422     0.831     1l|96 

				1|l|     1l|β_2   0.263     1.249     1l|72     0.168     0.971     1l|87     0.212     0.967     1l|80     0.170     0.966     1l|78     0.162     0.814     1l|79     0.181     0.710     1l|86     0.175     0.722     1l|87 

				1|l|     1l|β_3   2.953     1.092     1l|83     2.601     0.750     1l|91     3.207     0.756     1l|85     2.043     0.733     1l|89     1.998     0.691     1l|86     1.896     0.629     1l|93     1.599     0.780     1l|94                       
 1-24*Mental Health Score        1l|β_0        1.791     2.398     1l|86     1.659     2.092     1l|86      0.672     1.763     1l|88      0.356     1.616     1l|89      0.602     1.609     1l|85      0.367     1.614     1l|90      0.406     1.638     1l|87  

				1|l|     1l|β_1   2.101     3.778     1l|98     2.334     2.569     1l|100     1.178     1.587     1l|99      0.892     1.498     1l|98      0.952     1.487     1l|97      0.913     1.495     1l|98      0.968     1.515     1l|99  

				1|l|     1l|β_2   0.951     2.984     1l|82     1.042     2.492     1l|85      0.905     2.110     1l|84      0.651     2.004     1l|87      0.703     1.994     1l|88      0.659     2.012     1l|90      0.641     1.976     1l|89  

				1|l|     1l|β_3   1.338     3.619     1l|96     1.536     3.081     1l|100     1.459     2.164     1l|100     1.226     2.024     1l|100     1.285     2.054     1l|100     1.199     2.058     1l|100     1.272     1.988     1l|100   
 1-24*Physical Health Score        1l|β_0        0.695     2.158     1l|85     0.708     1.819     1l|85     0.343     1.579     1l|80     0.233     1.469     1l|86     0.309     1.452     1l|84     0.222     1.455     1l|88     0.197     1.488     1l|90 

				1|l|     1l|β_1   0.871     2.385     1l|95     0.898     1.294     1l|98     0.518     1.325     1l|97     0.457     1.247     1l|97     0.466     1.237     1l|97     0.425     1.230     1l|98     0.390     1.258     1l|98 

				1|l|     1l|β_2   0.684     2.640     1l|84     0.586     1.841     1l|88     0.378     1.823     1l|90     0.339     1.751     1l|92     0.391     1.742     1l|91     0.350     1.755     1l|93     0.331     1.731     1l|92 

				1|l|     1l|β_3   0.928     2.085     1l|87     1.139     2.079     1l|85     0.930     1.754     1l|85     0.810     1.656     1l|90     0.964     1.801     1l|87     0.876     1.654     1l|89     0.815     1.603     1l|92
 1-2
			tableA comparison of squared bias in point estimates of coefficients, standard deviations and 95% coverage of the five models of interest under seven treatments of missing data in the 100 sub sampled QASC data sets.
			
		



 §.§ Application to QASC Data Set

We now apply our proposed method to impute missing data in the original QASC data set with a total of 1480 patients. Unlike in Section 4.3 where we deliberately deleted some records so that we knew the true values, we do not know the true missing values here and therefore cannot measure imputation accuracy. We check the imputation quality by using diagnostics discussed in <cit.> and <cit.>.  Specifically, we examined the trace plots of the parameters and convergence in our proposed model (not shown here) and plotted the univariate densities/frequencies of the fully observed values (in black) and the average imputed values (in six colors) for some variables (see Figure 2). All the imputation methods generally agree with the complete data for the continuous variables `length of stay' and `age' and there are small disagreements for the variables `mental health score' and `physical health score'. The imputed values seem to be more spread out for `Bartell Index' than the observed data which is concentrated around 0. Overall, the frequencies of the categorical variables match the observed data with a few exceptions, for example, FCS imputes significantly more at level 4 for `Marital status'; and JM does not have any imputed values that fall into level 6 for `Modified Rankin Scale'. The departure from the observed data does not necessarily mean the imputation is poor, rather it may mean that the distribution of the missing data is different from what is observed, probably because of the missing data process is MAR rather than MCAR, lack of fit in the imputation model, etc. 




We also report the point estimates of coefficients as well as the standard deviations and p-values of the five models of interest in Table 7, by CC and the six imputation methods. While there are differences in the parameter estimates, the p-values across all the methods generally agree with each other, leading to the same clinical conclusions. There are some exceptions, for example in the random intercept logistic regression model for `Bartell Index 60', the coefficient of the interaction term β_3 is significant at the 0.1 level for the methods CC, FCS and our proposed Cluster Copula method, but significant only at the 0.05 level for the methods JM, Cluster JM, Cluster FCS and Copula_Hoff.

901
0.55

2|l|2*     3l|CC    3l|JM    3l|FCS    3l|Cluster JM     3l|Cluster FCS                                                            3l|Copula_Hoff    3l|Cluster Copula  
 3-23
2|l|                      point     sd     1l|p-value     point     sd     1l|p-value     point     sd     1l|p-value     point     sd     1l|p-value     point     sd     1l|p-value     point     sd     1l|p-value     point     sd     1l|p-value 
 4*Modified Rankin Scale 2        1l|β_0        -0.747                         0.156                       1l|.001                    -0.709                         0.152                       1l|.001                    -0.520                         0.147                       1l|.001                    -0.712                         0.160                       1l|.001                    -0.771                         0.150                       1l|.001                    -0.721                         0.151                       1l|.001                    -0.683                         0.154                       1l|.001                

  1|l|                           1l|β_1        0.114                          0.163                       1l|0.483                            0.091                          0.161                       1l|0.574                            0.179                          0.152                       1l|0.238                            0.102                          0.170                       1l|0.549                            0.131                          0.163                       1l|0.421                            0.092                          0.162                       1l|0.571                            0.125                          0.159                       1l|0.434                        

   1|l|                          1l|β_2        0.241                          0.226                       1l|0.285                            0.199                          0.218                       1l|0.362                            0.231                          0.213                       1l|0.279                            0.199                          0.224                       1l|0.375                            0.227                          0.216                       1l|0.293                            0.192                          0.221                       1l|0.384                            0.143                          0.228                       1l|0.530                        

    1|l|                         1l|β_3        -0.659                         0.244                       1l|0.007                            -0.588                         0.237                       1l|0.013                            -0.569                         0.225                       1l|0.011                            -0.601                         0.250                       1l|0.016                            -0.621                         0.242                       1l|0.010                            -0.616                         0.244                       0.012                            1l|-0.604                         0.240                       1l|0.012                        
1-24*Bartell Index 60        1l|β_0        2.349                          0.235                       1l|.001                    1.721                          0.181                       1l|.001                    2.436                          0.227                       1l|.001                    1.701                          0.167                       1l|.001                    1.829                          0.174                       1l|.001                    1.841                          0.182                       1l|.001                    2.010                          0.198                       1l|.001                

   1|l|                          1l|β_1        -0.200                         0.264                       1l|0.449                            -0.036                         0.216                       1l|0.869                            -0.402                         0.284                       1l|0.157                            -0.016                         0.209                       1l|0.939                            -0.134                         0.212                       1l|0.527                            -0.070                         0.228                       1l|0.759                            -0.058                         0.247                       1l|0.814                        

     1|l|                        1l|β_2        -0.472                         0.318                       1l|0.138                            -0.360                         0.248                       1l|0.146                            -0.494                         0.298                       1l|0.098                            -0.392                         0.236                       1l|0.096                            -0.466                         0.242                       1l|0.054                            -0.431                         0.245                       1l|0.078                            -0.369                         0.266                       1l|0.166                        

     1|l|                        1l|β_3        0.663                          0.380                      1l|0.081                            0.662                          0.311                       1l|0.033                            0.703                          0.375                       1l|0.061                            0.683                          0.301                       1l|0.023                            0.784                          0.327                       1l|0.016                            0.734                          0.323                       1l|0.023                            0.620                          0.337                       1l|0.066                        
 1-24*Bartell Index 95        1l|β_0        0.085                          0.161                       1l|0.598                            -0.098                         0.149                       1l|0.512                            -0.122                         0.149                       1l|0.413                            -0.125                         0.149                       1l|0.400                            -0.076                         0.140                       1l|0.591                            -0.065                         0.145                       1l|0.652                            -0.031                         0.149                       1l|0.838                        

   1|l|                          1l|β_1        0.288                          0.161                       1l|0.073                            0.258                          0.155                       1l|0.097                            0.216                          0.157                       1l|0.170                            0.282                          0.154                      1l|0.066                            0.241                          0.151                       1l|0.110                            0.312                          0.157                       1l|0.047                            0.305                          0.152                       1l|0.044                        

    1|l|                         1l|β_2        -0.163                         0.236                       1l|0.489                            -0.219                         0.214                       1l|0.304                            -0.197                         0.210                       1l|0.349                            -0.204                         0.217                       1l|0.346                            -0.211                         0.206                       1l|0.305                            -0.167                         0.213                       1l|0.434                            -0.144                         0.216                       1l|0.505                        

    1|l|                         1l|β_3        0.505                          0.242                       1l|0.037                            0.558                          0.232                       1l|0.016                            0.576                          0.228                       1l|0.012                            0.519                          0.228                       1l|0.023                            0.540                          0.230                       1l|0.019                            0.526                          0.231                       1l|0.023                            0.538                          0.230                       1l|0.020                        
 1-24*mental health score        1l|β_0        46.139                         0.810                       1l|.001                    45.106                         0.926                       1l|.001                    46.411                         1.206                       1l|.001                            44.656                         0.865                       1l|.001                    45.683                         0.833                       1l|.001                    45.174                         0.830                       1l|.001                    45.146                         0.813                       1l|.001                

1|l|                             1l|β_1        3.320                          0.931                       1l|.001                    3.426                          1.074                       1l|0.001                            2.476                          1.222                       1l|0.043                            3.826                          0.983                       1l|.001                    2.656                          0.963                       1l|0.006                            3.544                          0.977                       1l|.001                    3.259                          0.943                       1l|0.001                        

 1|l|                            1l|β_2        0.017                          1.201                       1l|0.989                           -0.404                         1.236                       1l|0.744                            0.002                          1.118                       1l|0.998                            -0.117                         1.212                       1l|0.923                            -0.509                         1.226                       1l|0.678                            -0.119                         1.255                       1l|0.924                            -0.235                         1.219                       1l|0.847                        

 1|l|                            1l|β_3        -0.067                          1.383                       1l|0.962                            0.671                          1.435                       1l|0.640                            -0.092                         1.386                       1l|0.947                            0.553                          1.385                       1l|0.690                            0.954                          1.403                       1l|0.497                            0.303                          1.426                       1l|0.832                            0.865                          1.339                       1l|0.518                        
 1-24*physical health score        1l|β_0        46.573                         0.840                       1l|.001                    45.145                         0.806                       1l|.001                    46.128                         0.951                       1l|.001                    45.035                         0.805                       1l|.001                    45.258                         0.786                       1l|.001                    45.383                         0.773                       1l|.001                    45.315                         0.802                       1l|.001                

   1|l|                              1l|β_1        -3.928                         0.808                       1l|.001                    -3.310                         0.846                       1l|.001                    -3.918                         0.892                       1l|.001                    -3.267                         0.833                       1l|.001                    -3.258                         0.836                       1l|.001                    -3.301                         0.803                       1l|.001                    -3.275                         0.819                       1l|.001                

     1|l|                            1l|β_2        -0.527                         1.226                       1l|0.667                             -0.481                         1.169                       1l|0.681                            -0.284                         1.128                       1l|0.801                            -0.518                         1.165                      1l|0.656                            -0.863                         1.131                       1l|0.446                            -0.645                         1.159                       1l|0.578                            -0.385                         1.155                       1l|0.739                        

    1|l|                             1l|β_3        3.035                         1.206                       1l|0.012                            3.148                          1.244                       1l|0.011                            2.944                          1.237                       1l|0.017                            3.343                          1.232                       1l|0.007                            3.505                          1.210                       1l|0.004                            3.340                          1.248                       1l|0.007                            3.015                          1.193                       1l|0.012                        
 1-2
tableA comparison of point estimates of coefficients, standard deviations and p-values of the five models of interest under seven treatments of missing data in the original QASC data set.







§ DISCUSSION

In this paper, we developed a copula based imputation model for multilevel data sets with mixed data. Copula based imputation models have a sound theoretical foundation and we have shown through simulations that copula based imputation models achieve reasonably accurate predictions of the missing values and recovery of parameters in some models of interest.

The copula based imputation models outperform the competing methods, especially when the variable distributions depart from normality. We also recommend taking into account clustering effects to incorporate information from the grouping structure in the analysis. This is confirmed from our simulation results, that when the ICC is high, imputation models with random effects added achieve better results.

One extension to our models is to add some `fixed' covariates. For the copula models in Section 3, all the variables appear on one side of the equations in (<ref>) and we model their relationship through the correlation matrices on the latent variable scale. But it is often of interest to see both the relationship among variables on the response side and the relationship between the responses and some covariates. For example in the QASC data set, `treatment' is fixed by design at hospital level, so we can treat it as a regressor. By doing so, the treatment effects on some process of care variables can be detected directly through the copula model on the latent variable scale. Here we consider variables with ordering, and extension to nominal variables is straightforward. Let i=1,...,m be the group index, j=1,...,n_i be the individual index within group i, and l=1,...,p be the variable index. Suppose the first k variables have common covariates x_i1,...,x_iq at the group level, in other words, they are fixed within group i. The correlation matrices for residual and random effects b_i are Γ and Ψ respectively as before, but the mean of the latent variables z is no longer zero. Again we use the monotone transformation z_ijl=Φ^-1(F(y_ijl)) to obtain the extended rank likelihood, then the model becomes:


    z_ij∼ N(b_i+x_i(β, 0), Γ), b_i∼ N(0,Ψ)
       ⇕
       (z_ij1,...,z_ijk,...,z_ijp)∼ N(b_i+(x_i1,...,x_iq)[ β_11    ⋯ β_1k    0    ⋯    0;    ⋮         ⋮    ⋮    ⋮; β_q1    ⋯ β_qk    0    ⋯    0 ],Γ))

It is straightforward to derive the full conditional distributions for the Gibbs sampler, and we omit the details here.

Choosing the form of copula is another issue which is a critical yet complicated task. <cit.>, <cit.> provide some guidance on choosing among existing copulas or creating new families of copulas. In this paper, we focused on the Gaussian copula because it is easy to extend to higher dimensions and computationally convenient. However, the main drawbacks of the Gaussian copula are the symmetry assumption and absence of tail dependence <cit.>. Therefore, some goodness-of-fit tests should be examined to check for a need to use other forms of copulas, for example, a (mixture of skewed) t-copulas.

kbib


42
	
	
	urlstyle[Abayomi et al.(2008)Abayomi, Gelman, and Levy]abayomi2008diagnostics
	Abayomi, Kobi; Gelman, Andrew, and Levy, Marc.
	Diagnostics for multivariate imputations.
	Journal of the Royal Statistical Society: Series C (Applied
		Statistics), 570 (3):0 273–291, 2008.
	
	[Aitchison and Bennett(1970)]aitchison1970polychotomous
	Aitchison, John and Bennett, Jo A.
	Polychotomous quantal response by maximum indicant.
	Biometrika, 570 (2):0 253–262, 1970.
	
	[Albert and Chib(1993)]albert1993bayesian
	Albert, James H and Chib, Siddhartha.
	Bayesian analysis of binary and polychotomous response data.
	Journal of the American Statistical Association, 880
	(422):0 669–679, 1993.
	
	[Bates et al.(2014)Bates, Mächler, Bolker, and
	Walker]bates2014fitting
	Bates, Douglas; Mächler, Martin; Bolker, Ben, and Walker, Steve.
	Fitting linear mixed-effects models using lme4.
	arXiv preprint arXiv:1406.5823, 2014.
	
	[Buuren and Groothuis-Oudshoorn(2011)]buuren2011mice
	Buuren, Stef and Groothuis-Oudshoorn, Karin.
	mice: Multivariate imputation by chained equations in r.
	Journal of Statistical Software, 450 (3), 2011.
	
	[Chib and Greenberg(1998)]chib1998analysis
	Chib, Siddhartha and Greenberg, Edward.
	Analysis of multivariate probit models.
	Biometrika, 850 (2):0 347–361, 1998.
	
	[Demarta and McNeil(2005)]demarta2005t
	Demarta, Stefano and McNeil, Alexander J.
	The t copula and related copulas.
	International Statistical Review/Revue Internationale de
		Statistique, pages 111–129, 2005.
	
	[Drechsler(2015)]drechsler2015multiple
	Drechsler, Jörg.
	Multiple imputation of multilevel missing data—rigor versus
	simplicity.
	Journal of Educational and Behavioral Statistics, 400
	(1):0 69–95, 2015.
	
	[Eddings and Marchenko(2011)]eddings2011accounting
	Eddings, W and Marchenko, Y.
	Accounting for clustering with mi impute.
	STATA. Available online at http://www. stata.
		com/support/faqs/statistics/clustering-and-mi-impute/, checked on,
	120 (5):0 2013, 2011.
	
	[Fox(2013)]fox2013package
	Fox, Maintainer John.
	Package ‘norm’.
	2013.
	
	[Genest et al.(1995)Genest, Ghoudi, and
	Rivest]genest1995semiparametric
	Genest, Christian; Ghoudi, Kilani, and Rivest, L-P.
	A semiparametric estimation procedure of dependence parameters in
	multivariate families of distributions.
	Biometrika, 820 (3):0 543–552, 1995.
	
	[Goldstein et al.(2009)Goldstein, Carpenter, Kenward, and
	Levin]goldstein2009multilevel
	Goldstein, Harvey; Carpenter, James; Kenward, Michael G, and Levin, Kate A.
	Multilevel models with multivariate mixed response types.
	Statistical Modelling, 90 (3):0 173–197,
	2009.
	
	[Graham et al.(2007)Graham, Olchowski, and Gilreath]graham2007many
	Graham, John W; Olchowski, Allison E, and Gilreath, Tamika D.
	How many imputations are really needed? some practical clarifications
	of multiple imputation theory.
	Prevention Science, 80 (3):0 206–213, 2007.
	
	[Gruhl et al.(2013)Gruhl, Erosheva, Crane,
	et al.]gruhl2013semiparametric
	Gruhl, Jonathan; Erosheva, Elena A; Crane, Paul K, and others, .
	A semiparametric approach to mixed outcome latent variable models:
	Estimating the association between cognition and regional brain volumes.
	The Annals of Applied Statistics, 70 (4):0
	2361–2383, 2013.
	
	[Hoff(2007)]hoff2007extending
	Hoff, Peter D.
	Extending the rank likelihood for semiparametric copula estimation.
	The Annals of Applied Statistics, pages 265–283, 2007.
	
	[Hollenbach et al.(2014)Hollenbach, Metternich, Minhas, and
	Ward]hollenbach2014fast
	Hollenbach, Florian M; Metternich, Nils W; Minhas, Shahryar, and Ward,
	Michael D.
	Fast & easy imputation of missing social science data.
	arXiv preprint arXiv:1411.0647, 2014.
	
	[Honaker et al.(2011)Honaker, King, Blackwell,
	et al.]honaker2011amelia
	Honaker, James; King, Gary; Blackwell, Matthew, and others, .
	Amelia ii: A program for missing data.
	Journal of Statistical Software, 450 (7):0
	1–47, 2011.
	
	[Käärik and Käärik(2009)]kaarik2009modeling
	Käärik, Ene and Käärik, Meelis.
	Modeling dropouts by conditional distribution, a copula-based
	approach.
	Journal of Statistical Planning and Inference, 1390
	(11):0 3830–3835, 2009.
	
	[Kole et al.(2007)Kole, Koedijk, and Verbeek]kole2007selecting
	Kole, Erik; Koedijk, Kees, and Verbeek, Marno.
	Selecting copulas for risk management.
	Journal of Banking and Finance, 310 (8):0
	2405–2423, 2007.
	
	[Kropko et al.(2013)Kropko, Goodrich, Gelman, and
	Hill]kropko2013multiple
	Kropko, Jonathan; Goodrich, Ben; Gelman, Andrew, and Hill, Jennifer.
	Multiple imputation for continuous and categorical data: Comparing
	joint and conditional approaches.
	Columbia University, Department of Statistics. New York, 2013.
	
	[Lascio(2015)]lascio2015exploring
	Simone GianneriniLascio, Alessandra Reale.
	Exploring copulas for the imputation of complex dependent data.
	Statistical Methods and Applications, pages 159–175, 2015.
	
	[Lee and Carlin(2010)]lee2010multiple
	Lee, Katherine J and Carlin, John B.
	Multiple imputation for missing data: fully conditional specification
	versus multivariate normal imputation.
	American Journal of Epidemiology, page kwp425, 2010.
	
	[Little and Rubin(2002)]little2014statistical
	Little, Roderick JA and Rubin, Donald B.
	Statistical analysis with missing data.
	John Wiley & Sons, 2002.
	
	[Liu and Wu(1999)]liu1999parameter
	Liu, Jun S and Wu, Ying Nian.
	Parameter expansion for data augmentation.
	Journal of the American Statistical Association, 940
	(448):0 1264–1274, 1999.
	
	[Lumley(2014)]lumley2014mitools
	Lumley, T.
	mitools: Tools for multiple imputation of missing data. r package
	version 2.0, 2014.
	
	[Middleton et al.(2011)Middleton, McElduff, Ward, Grimshaw, Dale,
	D'Este, Drury, Griffiths, Cheung, Quinn, et al.]middleton2011implementation
	Middleton, Sandy; McElduff, Patrick; Ward, Jeanette; Grimshaw, Jeremy M; Dale,
	Simeon; D'Este, Catherine; Drury, Peta; Griffiths, Rhonda; Cheung, N Wah;
	Quinn, Clare, and others, .
	Implementation of evidence-based treatment protocols to manage fever,
	hyperglycaemia, and swallowing dysfunction in acute stroke (QASC): a
	cluster randomised controlled trial.
	The Lancet, 3780 (9804):0 1699–1706, 2011.
	
	[Murray et al.(2013)Murray, Dunson, Carin, and
	Lucas]murray2013bayesian
	Murray, Jared S; Dunson, David B; Carin, Lawrence, and Lucas, Joseph E.
	Bayesian gaussian copula factor models for mixed data.
	Journal of the American Statistical Association, 1080
	(502):0 656–665, 2013.
	
	[Nelsen(2007)]nelsen2007introduction
	Nelsen, Roger B.
	An Introduction to Copulas.
	Springer Science & Business Media, 2007.
	
	[Pitt et al.(2006)Pitt, Chan, and Kohn]pitt2006efficient
	Pitt, Michael; Chan, David, and Kohn, Robert.
	Efficient bayesian inference for gaussian copula regression models.
	Biometrika, 930 (3):0 537–554, 2006.
	
	[Raghunathan et al.(2001)Raghunathan, Lepkowski, Van Hoewyk, and
	Solenberger]raghunathan2001multivariate
	Raghunathan, Trivellore E; Lepkowski, James M; Van Hoewyk, John, and
	Solenberger, Peter.
	A multivariate technique for multiply imputing missing values using a
	sequence of regression models.
	Survey Methodology, 270 (1):0 85–96, 2001.
	
	[Raghunathan et al.(2002)Raghunathan, Solenberger, and
	Van Hoewyk]raghunathan2002iveware
	Raghunathan, Trivellore E; Solenberger, Peter W, and Van Hoewyk, John.
	Iveware: Imputation and variance estimation software.
	Ann Arbor, MI: Survey Methodology Program, Survey Research
		Center, Institute for Social Research, University of Michigan, 2002.
	
	[Royston et al.(2005)]royston2005multiple
	Royston, Patrick and others, .
	Multiple imputation of missing values: update of ice.
	Stata Journal, 50 (4):0 527, 2005.
	
	[Rubin(1976)]rubin1976inference
	Rubin, Donald B.
	Inference and missing data.
	Biometrika, 630 (3):0 581–592, 1976.
	
	[Rubin(1987)]rubin1987multiple
	Rubin, Donald B.
	Multiple imputation for nonresponse in surveys (wiley series in
	probability and statistics).
	1987.
	
	[Schafer(1997)]schafer1997analysis
	Schafer, Joseph L.
	Analysis of Incomplete Multivariate Data.
	CRC press, 1997.
	
	[Schafer and Yucel(2002)]schafer2002computational
	Schafer, Joseph L and Yucel, Recai M.
	Computational strategies for multivariate linear mixed-effects models
	with missing values.
	Journal of computational and Graphical Statistics, 110
	(2):0 437–457, 2002.
	
	[Sklar(1959)]sklar1959fonctions
	Sklar, M.
	Fonctions de répartition à n dimensions et leurs
		marges.
	Université Paris 8, 1959.
	
	[Su et al.(2011)Su, Gelman, Hill, Yajima, et al.]su2011multiple
	Su, Yu-Sung; Gelman, Andrew; Hill, Jennifer; Yajima, Masanao, and others, .
	Multiple imputation with diagnostics (mi) in r: Opening windows into
	the black box.
	Journal of Statistical Software, 450 (2):0
	1–31, 2011.
	
	[Tanner and Wong(1987)]tanner1987calculation
	Tanner, Martin A and Wong, Wing Hung.
	The calculation of posterior distributions by data augmentation.
	Journal of the American Statistical Association, 820
	(398):0 528–540, 1987.
	
	[Trivedi and Zimmer(2007)]trivedi2007copula
	Trivedi, Pravin K and Zimmer, David M.
	Copula modeling: An Introduction for Practitioners.
	Now Publishers Inc, 2007.
	
	[Van Buuren(2007)]van2007multiple
	Van Buuren, Stef.
	Multiple imputation of discrete and continuous data by fully
	conditional specification.
	Statistical Methods in Medical Research, 160
	(3):0 219–242, 2007.
	
	[Zhao and Yucel(2009)]zhao2009performance
	Zhao, Enxu and Yucel, Recai M.
	Performance of sequential imputation method in multilevel
	applications.
	In American Statistical Association Proceedings of the Survey
		Research Methods Section, 2009.
	


