Introduction

In epidemiology, patient-related outcomes such as quality of life, wellbeing, or psychological health are usually measured using questionnaires. Mostly, these constructs encompass various subscales of the questionnaires. For instance the WHO-5 Well Being Index consists of subscales for positive mood, vitality, and general interest [1]. and the Depression Anxiety Stress Scales (DASS) contains a depression, an anxiety and a stress subscale [2]. Researchers typically report the following when using questionnaires with subscales: 1) the results of the analysis with the total score, and 2) the results of separate analyses performed for the different subscales. However, according to statistical theory this analytical approach is problematic, since subscales scores are correlated. In other words, when analysing the different subscales of a questionnaire in separate analyses, it is ignored that the subscales are filled in by the same subject and thus that subscale scores are not independent of each other. The widespread use of this problematic analytical approach is surprising, because multivariate multilevel analysis is a relatively simple alternative in which the dependency of observations is taken into account. In multivariate multilevel analysis, subscales of a particular questionnaire can be analysed in a single model, taking into account the correlation between subscales measured in the same subject. Although examples of medical studies using multivariate multilevel analysis for subscales exist [3-5], in most medical studies the dependency of the outcomes is ignored. Therefore, the present paper illustrates the use of multivariate multilevel analysis to analyse multiple correlated outcomes.

It should be noted that the use of multivariate multilevel analysis is not limited to questionnaires with different subscales. There are many more examples of studies in which multiple outcomes are measured in the same subjects (e.g., pathology in multiple brain regions [6], recovery of different segments of the heart after myocardial infarction [7], etc.). Because these outcomes are not independent of each other, also for analysing these type of data, multivariate multilevel analyses is highly suitable (see one of the examples used in this paper).

Multilevel analysis has been specifically developed to deal with correlated observations in educational research to take into account the dependency of the observations made on students who are in the same class (or even in the same school) [8]. Nowadays, multilevel analysis is used in many different research areas. For instance, in social epidemiology, multilevel analysis is used to take into account the dependency of the observations on subjects living in the same neighbourhood [9]. In clinical studies, multilevel analysis is, for instance, used for multicentre studies to taken into account the dependency of the observations on patients from the same centre [10,11]. Apart from suitability for crosssectional studies with correlated observations, multilevel analysis is also highly suitable for longitudinal studies in order to take into account the fact that the repeated observations within a subject are (highly) correlated [12]. Of note, the use of multilevel analysis for longitudinal data is actually the application of a multivariate multilevel analysis, as more than one outcome per subject is analysed in one model. In the literature, longitudinal data is always analysed as a multivariate outcome and separate analyses of the

different time-points are very rare. On the other hand, different outcomes from the same subject measured through for example subscales of a questionnaire are hardly ever analysed as a multivariate outcome. This is a rather peculiar phenomenon, because both longitudinal designs and designs with otherwise multiple outcomes per subject lead to highly correlated data, and this correlation should be taken into account in the analyses.

In sum, in the present paper the value of multivariate multilevel analyses for correlated outcomes was assessed by comparing multivariate multilevel analysis with separate analyses per outcome using both cross-sectional and longitudinal designs. Besides real life example datasets also simulated data was used in which various characteristics were systematically varied, such as the strength of the correlation between the multiple outcomes, the amount of missing data, and in longitudinal designs, the type of the covariate used in the analysis (i.e., time dependent versus time-independent).

Methods

In a multivariate multilevel analysis, a level for the different outcomes is present below the level of the subject (Figure 1). Repeatedly collecting multiple outcomes from one subject in a longitudinal study creates a three-level data structure. For this three-level structure basically two clustering options are possible: 1) The outcomes are clustered within the time-points and the time-points are clustered within the subjects or 2) The repeated measures are clustered within the outcomes and the outcomes are clustered within the subjects (Figure 2A & 2B).

Figure 1: A two-level multivariate multilevel model. Outcomes are clustered within subjects.

Figure 2A: A three-level multivariate longitudinal multilevel model. Outcomes are clustered within time-points, while the time-points are clustered within subjects.

Figure 2B: A three-level multivariate longitudinal multilevel model. Time-points are clustered within outcomes, while outcomes are clustered within subjects.

Cross-sectional real example dataset

The cross-sectional example dataset was taken from a single measurement round of a longitudinal observational cohort study, i.e., the Amsterdam Growth and Health Longitudinal Study (AGHLS) [13]. The aim of this longitudinal study was to follow up the natural growth, health, and lifestyle in a representative sample of 698 Dutch adolescents [13]. In total, ten measurement rounds were performed between 1976 and 2006. The cross-sectional example dataset was taken from the last measurement round in which subjective psychological well-being was measured with the 5-item World Health Organization Well-Being Index (WHO-5). The WHO-5 is focused on subjective quality of life based on three subscales: 1) positive mood, 2) vitality, and 3) general interest [1]. In the example, the difference between males and females for the three different subscales of wellbeing was investigated. The multivariate cross-sectional multilevel model contained two-levels, i.e., the three subscales of well-being were clustered within the subjects (Figure 1).

Longitudinal real example datasets

The first longitudinal example dataset was taken from the Netherlands Obsessive Compulsive Disorder Association (NOCDA) study. NOCDA is a longitudinal cohort study investigating the long-term course of obsessive-compulsive disorder (OCD) in patients referred to mental health care centres [14]. In this study, 419 patients were followed for a period of 6 years. Baseline measurements were performed between 2005 and 2009 and after that, three follow-up measurements were performed after 2, 4 and 6 years. In the present example, data from the three follow-up measurements were used on two related outcomes: depression and anxiety. The second longitudinal example dataset is taken from the last three measurement rounds of the AGHLS (see above) in which the correlated outcomes were three risk factors for cardiovascular disease (i.e., triglycerides, body mass index (BMI), and mean arterial blood pressure). Because the risk factors were measured on totally different scales, z-scores were used for the analyses (z-scores were taken for each outcome separately, across measurements, to allow for detection of both between subject and within subject [i.e., time] effects).

In both longitudinal example datasets, three analyses were performed. 1) the development over time, 2) the longitudinal relationship with a time-independent covariate (i.e., sex in both examples), and 3) the longitudinal relationship with a timedependent covariate (i.e., symptom severity in the first example and the personality characteristic dominance in the second example). For the analyses, two multivariate longitudinal multilevel models were used, both containing three levels: 1) the outcomes were clustered within the three time-points, while the three time-points were clustered within the subjects, 2) the repeated measures were clustered within the outcome and the outcomes were clustered within the subject (see Figure 2a and 2b respectively). In all examples, the results of the multivariate multilevel analyses were compared to the results of the analyses for the outcomes separately.

All analyses were performed with STATA (version 17).

Simulations

For each simulation scenario, 1000 samples were generated with a continuous normally distributed outcome variable. In the cross-sectional simulated dataset, three continuous outcome variables mimicked a questionnaire with three subscales. Furthermore, a dichotomous covariate was created equally divided over the subjects. In the cross-sectional analysis the relationship between the dichotomous covariate and the different subscales of the outcome questionnaire was analysed. This was done with 1) a multivariate cross-sectional multilevel analysis with two levels: subscales were clustered within the subject and 2) three linear regression analyses for the three subscales separately. Different scenarios were created: different standard deviations for the three subscales, different correlations between the three subscales and different number of subjects, which was either 100 or 200 (resulting in a total number of 300 and 600 observations, respectively). Lastly, several datasets with increasing amount of missing data (missing data percentages ranging between 10% and 40%) and different types of missing data (missing completely at random [MCAR] and missing at random [MAR]) in some of the subscales were simulated.

The longitudinal simulated dataset consisted of 100 subjects who were measured at three time-points. In this dataset, the continuous outcome variable were z-scores of three risk factors (i.e., comparable to the second longitudinal real example dataset, so z-scores were created for each outcome separately, across measurements). This led to a total number of 900 observations. Three analyses were performed: 1) the development over time, 2) the relationship with a time-independent covariate and 3) the relationship with a time-dependent covariate. All three were analysed with 1) a longitudinal multivariate multilevel analysis with three levels: outcomes clustered within the time-point and subsequently the repeated measures clustered within the subject 2) a longitudinal multivariate multilevel analysis with three levels: repeated measures clustered within the outcome and subsequently the outcomes clustered within subjects and 3) three separate longitudinal multilevel analyses with two levels: the three repeated measures clustered within the subject. The development over time and the relationships with the time-dependent and timeindependent covariates were simulated to be different for the three outcomes. In the longitudinal simulations, different scenarios were created with different correlations between the outcomes and different correlations within an outcome over time (Table 1). Furthermore, several datasets with increasing amount of missing data at the repeated measures (i.e., 25% missing at the second measurement and 50% missing at the third measurement) and with different types of missing data (i.e., MCAR and MAR)) were simulated. All simulations were performed with STATA (version 17). For both the cross-sectional and longitudinal simulations for MCAR, missing data was randomly drawn from the whole population, while for MAR, missing data was randomly drawn from the lower 75% of the distribution of the outcome. It should be noted that the aim of the simulation study is only to compare the different methods with each other under changing circumstances.

Table 1: Different correlations used for the longitudinal simulations.

(Tables 2 and 3) show the data structure needed for both the cross-sectional and the longitudinal multivariate multilevel analysis. (Box 1 and Box 2) show examples of the STATA code for both the cross-sectional and longitudinal simulations.

Table 2: Data structure used for the cross-sectional multivariate multilevel analysis.

Table 3: Data structure used for the longitudinal multivariate multilevel analysis.

Box 1: Example of the STATA code for one of the cross-sectional simulations.

Box 2: Example of the STATA code for one of the longitudinal simulations.

Results

Cross-sectional example dataset

(Table 4) shows descriptive information regarding the cross-sectional example dataset regarding subjective psychological well-being in the 2006 measurement of the AGHLS study.

Table 4:  Mean and standard deviation ( ) of the three subscales of wellbeing for males and females

(Table 5) shows the results of the analyses related to the differences between males and females in three subscales of wellbeing. It should be noted that within multilevel analysis, the dependency of the observations is usually quantified by the intraclass correlation coefficient (ICC) [15]. In the cross-sectional multivariate multilevel analysis, there is one ICC (i.e.,on subject level), which reflects the average correlation between the outcomes within a subject.

Table 5: Results* of the analyses related to the difference between females and males in three subscales of wellbeing; ICC = intraclass correlation coefficient, *regression coefficients and standard errors ( )

The results show that the regression coefficients indicating the difference between males and females as well as the standard errors of the regression coefficients were exactly the same for the three separate regression analyses and the multivariate multilevel analysis. This, even though, the ICC on subject level was relatively high.

Longitudinal example datasets

(Table 6) shows descriptive information of the longitudinal example datasets, regarding the development of depression and anxiety across three measurements in the NOCDA study (Example 1) and the development of three risk factors for cardiovascular disease (triglycerides, BMI and blood pressure) across three measurements in the AGHLS study (Example 2). In a longitudinal multivariate dataset, there are two ICCs. One reflecting the average correlation between the outcomes within the subject per time-point (i.e., the ICC over outcomes), and one reflecting the average correlation between the measurements over time (i.e., the ICC over time).

Table 6:  Mean and standard deviation ( ) of the outcomes of the two longitudinal example datasets at the three follow-up measurements; ICC = intraclass correlation coefficient; *ICC based on z-scores across measurements.

From Table 6 it can be seen that in both datasets the number of subjects decreased over time, which is typically seen in longitudinal studies. Regarding the development over time, in the first example dataset there was a sharp decrease from measurement 1 to measurement 2 in both outcomes, while in the second dataset there were only relatively small changes over time. It can also be seen that in the first longitudinal example dataset the ICCs over outcomes within a time-point were slightly higher than the ICCs over time, while in the second longitudinal example dataset, the ICCs over outcomes were much lower than the ICCs over time.

Before performing the various analyses, it should be evaluated which of the two longitudinal multivariate multilevel models should be used in the examples by assessing which model best fits the data. Therefore, for both examples, an intercept only model was analysed. (Table 7) shows the ICCs and the Bayesian Information Criteria (BIC) obtained from the different analyses on the two longitudinal example datasets. The BIC is a model fit indicator which is often used to compare models with each other. A lower BIC value indicates a better model fit [16]. Regarding the first longitudinal multivariate multilevel analysis, the ICC on time level refers to the correlation between the outcomes at a particular time-point, while the ICC on subject level refers to the correlation between the repeated measures of an outcome within a subject. Regarding the second longitudinal multivariate multilevel analysis, the ICC on outcome level refers to the correlation between the repeated measures within an outcome, while the ICC on subject level refers to the correlation between the outcomes within the subject. From the BIC values, it is obvious that for example 1, the first longitudinal multivariate multilevel model fitted best, while for example 2, the second longitudinal multivariate multilevel model fitted best. In the remaining part of the analyses on the example datasets, these best fitting models were used (Table 8).

Table 7: Intraclass correlation coefficients and model fit of the different multivariate multilevel analysis.

Table 8: Results* of the analyses related to the development over time in different outcomes;*regression coefficients and standard errors ( )

The results of the different analyses investigating the development over time in the two outcomes (depression and anxiety) of the first longitudinal example dataset show that there was a small difference in both the estimated regression coefficients, and the standard errors between the longitudinal multivariate multilevel model and the two separate longitudinal multilevel models. For the second longitudinal example dataset, again small differences in estimated regression coefficients were found. The standard errors obtained from the longitudinal multivariate multilevel analysis on the other hand were lower than the one obtained from the separate analyses for triglycerides and mean blood pressure, but higher than the one obtained from the separate analysis for BMI.

(Table 9) shows the results of the different analyses investigating the difference in outcomes between males and females on average over time, and (Table 10) shows the results of the different analyses investigating the relationship between symptom severity (in the first longitudinal example dataset) and dominance (in the second longitudinal example dataset) and the outcomes on average over time.

Table 9: Results* of the analyses related to the difference between females and males in the different outcomes on average over time;