{'fraction_non_alphanumeric': 0.024609193623278132, 'fraction_numerical': 0.01501315585822628, 'mean_word_length': 4.824368990384615, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Employing latent variable models to improve efficiency in composite endpoint analysis

Composite endpoints that combine multiple outcomes on different scales are common in clinical trials, particularly in chronic conditions. In many of these cases, patients will have to cross a predefined responder threshold in each of the outcomes to be classed as a responder overall. One instance of this occurs in systemic lupus erythematosus (SLE), where the responder endpoint combines two continuous, one ordinal and one binary measure. The overall binary responder endpoint is typically analysed using logistic regression, resulting in a substantial loss of information. We propose a latent variable model for the SLE endpoint, which assumes that the discrete outcomes are manifestations of latent continuous measures and can proceed to jointly model the components of the composite. We perform a simulation study and find the method to offer large efficiency gains over the standard analysis. We find that the magnitude of the precision gains are highly dependent on which components are driving response. Bias is introduced when joint normality assumptions are not satisfied, which we correct for using a bootstrap procedure. The method is applied to the Phase IIb MUSE trial in patients with moderate to severe SLE. We show that it estimates the treatment effect 2.5 times more precisely, offering a 60% reduction in required sample size.


Introduction
Composite endpoints combine multiple outcomes in order to determine the effectiveness or efficacy of a treatment for a given disease. They are typically recommended when a disease is complex or multi-system and meaningful improvement cannot be captured in a single outcome.
Furthermore, the endpoint may be a combination of continuous and discrete outcomes which are collapsed in to a single binary responder index. Table 1 shows examples of diseases that use composite endpoints combining multiple continuous and discrete components. Responders in fibromyalgia must respond in two continuous and one ordinal component however responders in trials for frailty or soft tissue infections must respond in a total of five continuous and discrete components. Generally, these composite responder endpoints will be treated as a single binary outcome and analysed using a logistic regression model, which we term the standard binary method. This solves problems with multiplicity however results in large losses in efficiency (Wason and Seaman (2013)). The aim of this paper is to propose a joint modelling framework within which we can model the components of the composite, retaining the information on the original scales of the outcomes, hence increasing efficiency. One likelihood based method for handling mixed data is the factorisation model. The objective is to factorise the  (2013)). This accounts for correlations between the outcomes by including one response as a covariate in the model for the other response. In the graphical modelling literature this has been termed the 'Conditional Gaussian Distribution' (Whittaker (1990); Lauritzen and Wermuth (1989)). An advantage of these methods in relation to the composite endpoint problem is that we may account for correlations between measurements whilst making inference directly on the outcomes that we have measured, hence they fall within a broader class of 'direct methods'. Examples of applications of these ideas, which build on the work of Olkin and Tate (1961), include developmental toxicity studies by Fitzmaurice and Laird (1995) and in the longitudinal setting, the augmented binary method was developed for application to composite endpoints in clinical trials where the composite is formed of one continuous and one binary outcome (Wason and Seaman (2013); Wason and Jenkins (2016); McMenamin and others (2018)). One difficulty with these methods beyond the bivariate scenario is the range of possibilities for the factorisations, with no consensus on how this should be determined. In the case of the SLE responder endpoint with four components, this amounts to 24 possible factorisations, each of which may result in different conclusions (Verbeke and others (2014); de Leon and Carriere (2013)).
Another family of models used to model mixed outcome types which feature frequently in economics and finance are copulas. These are functions that join or couple multivariate distribution functions to their uniform one-dimensional marginal distribution functions as discussed by Nelsen (1999). Copulas offer a flexible framework in this setting, as the marginal distribution functions need not come from the same parametric family. While the construction of copulas is considered to be mathematically elegant and the flexibility with which we can model appealing, they are not without their shortcomings. Extensions beyond the bivariate setting are difficult and have failed to perform well in many applications (de Leon and Carriere (2013)). Other practical implications include poor out-of-sample predictions due to the wide variety of copulas available. These restrictions, along with difficulties in longitudinal settings with unbalanced data structures, have seen few applications of copulas for mixed outcome types in the medical statistics literature (Verbeke and others (2014)). Applications of copulas in mixed outcome settings include de Leon and Wu (2010) and Wu and de Leon. A.R. (2014).
Another likelihood based method that allows for more flexibility when modelling the correlations between outcomes falls within the framework of latent variable models (Skrondal and Rabe-Hesketh (2004)). The multiple outcomes are assumed to be physical manifestations of some underlying latent process. This is modelled by including the same latent variable in each of the models for the observed responses. The outcomes are then assumed to be independent conditional on this latent variable. This solves the problem of deciding the order of factorisations in previously discussed methods however this formulation results in the inclusion of some covariance parameters in the mean structure, leaving the model sensitive to misspecification of the correlation structure (Sammel and Ryan (2002)). One example of these models is introduced by Sammel and others (1997), where effects of covariates of interest are modelled through this shared latent variable. Although these models have the intuitive interpretation that each outcome is attempting to capture underlying disease activity, the correlation matrix is restricted to allow for the same correlation between each pair of outcomes, which is unlikely in practice. This structure is relaxed in work by Dunson (2000), where the effects of covariates are included in the model separately from the latent variable. The correlation structure can be further relaxed to allow for a different latent variable for each outcome, meaning that pairs of outcomes are not assumed to have the same correlation. However these models would require integrating out each of the latent variables in order to obtain the joint distribution of interest (McCulloch (2008)). Furthermore, they are relevant in applications with multiple time points however less so for a single time point, as is the case for the composite endpoint problem.
Latent variables have also been used in the setting of mixed continuous and discrete variables to a different end. Namely, the outcomes adopt a correlated Gaussian distribution by assuming that the discrete outcomes are coarsely measured manifestations of underlying continuous variables subject to some threshold specifications, as seen in Ashford and Sowden (1970) and Chib and Greenberg (1998). Specifying discrete variables in terms of a partitioning of the latent variable space into non-overlapping intervals dates back to Pearson (1904) in relation to his generalised theory of alternative inheritance and has received much consideration in the literature since. Terminology surrounding these models is inconsistent but they are often referred to as multivariate probit models (Ashford and Sowden (1970)). This theory can also be found to underpin conditional grouped continuous models (Poon and Lee (1987)). By formulating the distribution in this way, we can correlate the error terms between models and work within the familiar paradigm of Gaussian distributions and maximum likelihood theory. The theory and application of these ideas for a mixture of continuous and binary outcomes has featured in the statistics literature, see for example work by Tate (1955), Cox and Wermuth (1992) and Catalano and Ryan (1992).
Generalisations of these ideas, which appear less frequently in the literature, lead to methods for modelling continuous and ordinal variables, with applications in developmental toxicology (Catalano (1997); Samani and Ganjali (2008); Arminger and Kusters (1988); Regan and Catalano (2000); Faes and others (2002)). Despite the advantages, the multivariate probit model has not realised its full potential in the applied biostatistics literature. This was noted by Lessafre and Molenberghs (1991) and we believe it still to be the case today. The few applications that do appear tend to demonstrate bivariate scenarios or those that mix continuous and binary or continuous and ordinal, rather than all three. Furthermore, these models have not been considered specifically to address challenges with modelling composite endpoints. Other work has combined thresholding the response variables and introducing latent variables in the model, examples include work by Gueorguieva and Agresti (2001), de Leon and Carriere (2013), Gueorguieva and Sanacora (2003) and Gueorguieva and Sanacora (2006), however these ideas are most applicable in the longitudinal setting. We will therefore employ the latent outcome framework and investigate its use for the composite endpoint problem.
The paper proceeds as follows. In Section 2 we discuss systemic lupus erythematosus (SLE), the motivating example for the methods. In Section 3 we introduce the latent variable model and discuss how we conduct estimation and inference. In Section 4 we compare the behaviour of the latent variable model with the augmented binary and standard binary methods, including the case when the key assumptions are not satisfied. In Section 5 we apply the methods to the Phase IIb MUSE trial in patients with moderate to severe SLE. Finally, in Section 6 we discuss our findings and make some recommendations.

Motivating example
In what follows we will focus specifically on systemic lupus erythematosus (SLE), however the methods introduced will be relevant to other diseases using endpoints with a similar structure. The SLE endpoint is shown in Figure 1. It combines a continuous PGA measure, a continuous SLEDAI measure, an ordinal BILAG measure and a binary corticosteroids measure, where patients must meet the response criteria in all components in order to be classed as a responder overall. Note that the SLEDAI and BILAG measures are themselves composite scores deriving from a combination of items, however this will not be considered in the analysis.
The real data underpinning this motivation comes from the MUSE study (Furie and others (2017)). It was a Phase IIb, randomised, double-blind, placebo-controlled study investigating the efficacy and safety of anifrolumab in adults with moderate to severe SLE. Patients (n=305) were randomised to receive anifrolumab (300mg or 1000mg) or placebo, in addition to standard therapy every 4 weeks for 48 weeks. The primary end point was the percentage of patients achieving an SRI response at week 24 with sustained reduction of oral corticosteroids (<10mg/day and less than or equal to the dose at week 1 from week 12 through 24). The methods discussed will make inference at one time point, as this is the case in the trial, although they can be easily extended for the longitudinal case.

Notation
. T i represents the treatment indicator for patient i, y i10 and y i20 are the baseline measures for Y i1 and Y i2 respectively.

Model
The mean structure for the outcomes is shown in (3.1). The baseline measures y 10 and y 20 are included in the model for Y 1 and Y 2 respectively.
The observed discrete variables are related to the latent continuous variables by partitioning the latent variable space, as shown in (3.2). The lower and upper thresholds for both discrete variables are set at τ 03 = τ 04 = −∞, τ 53 = τ 24 = ∞. The intercept term for the ordinal variable in (3.1) is set at γ 0 = 0 so that the cut-points τ 13 , τ 23 , τ 33 , τ 43 may be estimated. The intercept for the binary outcome ψ 0 may be estimated, as τ 14 = 0.
Following these assumptions, we can model the error terms in (3.1) as multivariate normal with zero mean and variance-covariance matrix Σ, as shown in (3.3). Note that the error variances for ε * 3 , ε * 4 are σ 3 = 1 and σ 4 = 1. This does not represent a constraint on the model but rather a rescaling required for identifiability.
where θ is a vector which contains all model parameters. The observed likelihood can then be expressed as in (3.5). I{Yi3=w,Yi4=k} (3.5) The joint probability of patients having discrete measurements Y i3 = w and Y i4 = k must be multiplied over the five ordinal levels and two binary levels resulting in ten combinations of the probabilities in (3.6) to be calculated. We discuss the intuition for (3.6) in Appendix A.

Estimation
As the variance parameters (σ 1 , σ 2 ) are required to be greater than 0, we introduce parameters (δ 1 , δ 2 ) such that σ 1 = exp(δ 1 ) and σ 2 = exp(δ 2 ). This transformation ensures that the variance is above 0 whilst allowing the parameter we estimate to take any real value. We must also ensure that the correlation parameters (ρ 12 , ρ 13 , ρ 14 , ρ 23 , ρ 24 , ρ 34 ) are estimated within (-1,1) by introducing (δ 12 , δ 13 , δ 14 , δ 23 , δ 24 , δ 34 ), where ρ 12 = 2expit(δ 12 ) − 1, We fit the model in R by coding the likelihood function, probability of response and using the delta method to obtain standard errors. The bivariate distribution functions in (3.6) are estimated using 'pmvnorm', using the method of Genz (1992). The likelihood maximisation is conducted using the 'nlminb' function in the 'optimx' package, which is the best performing method in terms of accuracy and convergence rate, however is the slowest. We use the 'Hessian' function in the 'numDeriv' package to obtain the Hessian matrix and invert this to get the covariance matrix of the model parameters. In a small number of cases the Hessian is not positive definite because of computational error, meaning that it cannot be inverted. This is rectified in these cases by using the 'near PD' function in the 'Matrix' package, which computes the nearest positive definite matrix.

Inference
We wish to make inference on the probability of response. Let S i be an indicator for patient i denoting whether or not they achieved response defined by We obtain the integrand in (3.8) by using the fitted values of the parameters in the conditional mean and conditional covariance matrix in (3.7). Parameter estimates from these methods are maximum likelihood estimates and so we avail of asymptotic maximum likelihood theory. The integral in (3.8) is evaluated using the 'R2Cuba' package to obtain estimates for each patient, assuming they were treatedp i1 and not treatedp i0 . The odds ratio treatment effect is then defined as shown in (3.9).δ Note that we can easily define a risk difference or risk ratio using these quantities but in what follows we considerδ to be the effect of interest. The standard error estimates are obtained using the delta method. This requires the covariance matrix of the maximum likelihood estimates Cov(θ) and δ , the vector of partial derivatives ofδ with respect to each of the parameter estimates. The variance ofδ is obtained as shown in (3.10).
Alternatively, the quantity in (3.8) can also be considered to be a multivariate Gaussian hidden truncation distribution, from which we can obtain a closed form solution, and proceed as detailed by Arnold (2009).
Another important consideration for the model is how to assess goodness-of-fit. We propose an extension to an existing method for application in this case, which is detailed in Appendix B in the supplementary material.

Simulation study
We are interested in comparing the performance of the latent variable, augmented binary and standard binary methods through simulation. The models for the augmented binary and standard binary methods are included in Appendix C in the supplementary material.
The methods are evaluated against a range of performance criteria, which are included with their Monte Carlo standard errors in Appendix E of the supplementary material. For further details see Morris and others (2017). Figure 2 shows the bias of the methods as the treatment effect varies. The standard binary method is unbiased, as we would expect for a logistic regression in a large sample. The latent variable method is unbiased for smaller treatment effects but a small bias towards the null is introduced as the treatment effect increases. The augmented binary method is biased away from the null in this setting and the bias increases as the treatment effect increases. Given that this performance is worse than is suggested from previous applications of the augmented binary method in Wason and Seaman (2013) and Wason and Jenkins (2016), this would suggest that the treatment effect from the augmented binary method may be biased if the model is misspecified. The coverage of the methods is shown in Figure 3. The binary method has  Fig. 3. Coverage probability reported from the latent variable method, augmented binary method and standard binary method for nsim=5000, total sample size N=300 for true log-odds treatment effect between 1.2 and 2.2. The composite endpoint of interest contains four components: two continuous, one ordinal, one binary and treatment effects are present in all four components we can look at bias-corrected coverage, as recommended in Morris and others (2017). Figure 4 shows both the coverage and bias-corrected coverage for the three methods. The properties of the standard binary method remain unchanged. The bias-corrected coverage of the latent variable method is 0.95, which indicates that any under-coverage is due to the bias present. This is not true for the augmented binary method which shows small improvements in bias-corrected coverage, indicating that under-coverage is present in this method due to reasons other than bias. Again, this may be down to model misspecification. The power of the three methods is shown in Figure   5. The performance of the binary and augmented binary methods are as we would expect based on previous findings in Wason and Seaman (2013) and Wason and Jenkins (2016 that the standard binary method has the smallest bias and that the latent variable method has the smallest variance. The mean squared error (MSE) provides a combined measure of bias and variance. Figure 6 shows the MSE of the three methods as the treatment effect varies. The MSE for the standard and augmented binary methods is approximately 6.5 times that of the latent variable method. However, this measure should be interpreted with care due to the fact that the MSE is more sensitive to the sample size than comparisons of bias or empirical SE alone (Morris and others (2017)). To understand more about the precision performance of the augmented binary method in particular, we vary the responder threshold θ 1 to change the proportion of responders in that outcome. Figure 7 shows the density of the Y 1 variable and the relative precision of the methods, as the responder threshold varies. The precision gains from the augmented binary method diminish as the threshold increases. This is intuitive, as improvements in efficiency fall as the continuous component becomes less responsible for driving response. It is interesting to note that all precision gains are lost for any thresholds above -4. Therefore, even when 20% of patients are non-responders, all efficiency gains are lost. The percentage of responders needed to improve efficiency using the augmented binary method will of course depend on the correlation structure employed. Due to the additional information in the other components, the latent variable method is still five times as precise as the other methods. driving response. Figure 8 shows boxplots of the relative precision for the methods for four different response combinations, namely when response is driven by

Varying treatment effect
and (Y 4 ), where Y 1 and Y 2 are observed as continuous variables, Y 3 is ordinal and Y 4 is binary.
When all four components contribute to response, the latent variable method outperforms the other methods, offering large precision gains. The variability in the magnitude of these gains is large, with the median result showing that the latent variable method reports the treatment effect 8 times more precisely than the binary method and 6 times more precisely than the augmented binary method. If response is driven by (Y 1 , Y 2 , Y 3 ) then the relative median gains for the latent variable method are larger, however note that in less that 2% of cases the treatment effect is reported equally or less precisely than from both of the other methods. The findings are similar when response is driven by (Y 1 , Y 4 ), however the median gains are much smaller. The treatment effect is reported 5 times more precisely from the latent variable method than the binary method in this setting. Note that as the augmented binary method models the relevant components it still performs well and again better than the latent variable method in a very small number of cases.
When binary Y 4 determines response, the augmented binary method offers no improvement in precision whereas the latent variable method is approximately 1.5 times more precise. It is clear from the results that the magnitude of the precision gains from the latent variable method is highly dependent on the structure of the data.

Sensitivity analysis
The key assumptions in this model are that of joint normality of the four components and that the discrete variables are realisations of latent continuous variables. Although it is not possible to test these assumptions in real data, we can investigate how robust the latent variable method is to deviations from these conditions. We can do this by drawing from the multivariate skew-

Data structure
Due to data sharing policy, we conduct the analysis for a subset of the patients, N=278 rather than N=305 reported in the paper, so the results will differ from the original paper. Furthermore, only the anifrolumab 300mg arm (n=95) and the placebo arm (n=87) will be used to illustrate the methods.
The simulation results have suggested that the structure of the data is important for how the methods will perform, in particular the magnitude of the precision gains depends highly on which components drive response. Table 3 shows the criteria for response in each component and the rates of response in each by treatment arm. This suggests that the components responsible for responder discrimination are the continuous SLEDAI measurement and the binary taper measure.
We can further explore the structure of the data by visualising the 4-D endpoint. Figure 9 shows a plot of the four components in the SLE index. The two panels show taper responders and non-responders, the levels in BILAG are denoted using colours where any coloured data points representing Grade B -Grade E are responders. The response thresholds for the continuous measurements are included, where a patient must be below the threshold to be considered a responder. We can conclude that response is entirely driven by SLEDAI and the taper variable, as there are no PGA non-responders not already accounted for by SLEDAI and no purple data points in the responder quadrant.  SLEDAI is plotted on the y-axis and PGA on the x-axis, along with their corresponding dichotomisation thresholds. Levels of BILAG are represented by different colours and taper responders and nonresponders are split across two panels.

Results
The probability of response in the placebo arm is estimated as 0.199 by the latent variable method, 0.211 by the augmented binary method and 0.224 by the standard binary method.
A much larger discrepancy between the methods is shown in the treatment arm, where the probability of response is estimated at 0.311, 0.324 and 0.382 in the latent variable, augmented binary and standard binary methods respectively.
The log-odds treatment effect point estimates and confidence intervals for the MUSE trial are shown in Table 4. Both joint modelling methods estimate the treatment effect more precisely.
Although there may be bias present in the point estimates for the joint modelling methods, the confidence intervals entirely overlap with that of the binary method. All three methods indicate that anifrolumab 300mg performs better than placebo, as in the original findings. The latent variable model fits the data well according to the modified Pearson residuals, see Appendix G.
The simulation results indicated that the latent variable method may report the treatment effect with bias and have problems with bias related under-coverage when the treatment effect is large and when the assumption of joint normality is not satisfied. As the problems with the performance are bias related, we suggest implementing a bootstrap procedure to correct for this.
In this scenario N=182 and n boot =1000, therefore the procedure is as follows:  In terms of estimated precision, it is interesting to determine where the trial data set lies in the distribution of datasets generated in the simulation study. The latent variable method reports the treatment effect 2.5 times more precisely than the standard binary method in this setting, whilst the augmented binary method is 2.4 times more precise. We would have expected this similar performance as the augmented binary method models the SLEDAI and taper variablesthe only components driving response. This increase in precision from the latent variable method compared with the binary method amounts to a 60% reduction in required sample size.

Discussion
In this paper we addressed the issue of substantial losses of information when modelling complex composite endpoints. By employing concepts of partitioning latent variable outcome spaces we could model the observed structure of the composite endpoint, which resulted in large gains in efficiency. Sensitivity analyses showed that a bias is introduced when the assumptions of joint normality were not satisfied, however similar reductions in variance were observed. When applying the methods to the MUSE trial, we implemented a bootstrap procedure to correct for the presumed bias, as joint normality could not be assessed. The treatment effect was reported 2.5 times more precisely than that reported from the standard binary method.
Bias correction appears to perform well in the real data, where the crucial assumptions cannot be tested. The point estimate is shifted by a magnitude that would have been expected from the simulation results and the estimate of the variance is similar to that obtained in the single trial dataset. Furthermore the bootstrap confidence interval for the treatment effect is contained within that for the binary method, which offers further reassurance for application. However, more work could be done to investigate different structures and scenarios to ensure that the bias correction is always performing as expected. Ideally, we would investigate this further across a large number of datasets however this is too computationally intensive. To perform this on one replicate, where n boot = 1000 using 200 cores on a high performance computer (HPC) currently takes 7 hours. Exploring this further through bootstrapping or employing alternative multivariate distributions is an area for future research.
The precision gains offered by the latent variable method offer justification for the additional complexity. However, the magnitude of these gains are highly dependent on the components that drive response. The baseline case in the simulations was chosen to reflect when a composite endpoint is recommended for use, i.e. when all four components determine response rates. In this scenario, the precision gains achieved resulted in the latent variable method reporting the effect 2.5 to 17.5 times more precisely than the standard binary method. However, in practice in SLE trials, this has not been found to be the case. A review of two phase 3 trials (N= 2262) using the SRI-5 index found the SRI-5 response rate at week 52 for all patients was 32.8% (Kalunian and others (2018)). Non-response due to a lack of SLEDAI improvement, concomitant medication noncompliance or dropout was 31, 16.5 and 19.1%, respectively. Non-response due to deterioration in BILAG or Physicians Global Assessment after SLEDAI improvement, concomitant medication compliance and trial completion was 0.5%. This is in agreement with our findings from the MUSE trial data, which suggests that the precision gains in the baseline case are optimistic. The simulation results show that when one continuous and one binary component drive response, the latent variable method may be anywhere between 1 and 12 times as precise as the binary method and up to 7 times as precise as the augmented binary method. In a very small number of cases (<2%) there are no efficiency gains from using the latent variable method in this scenario. However the potential gains available in 98% of cases ensures that implementing the latent variable method is still very much a worthwhile endeavour, for all stakeholders in a clinical trial.
In addition to SLE, we have identified other disease areas that have a similar complex composite structure, meaning the potential to improve efficiency extends well beyond SLE. However, it must be acknowledged that the exact structure of the endpoint may offer different magnitudes of bias and precision, and may require longer computational time. Furthermore, in conditions where longitudinal data is required to sufficiently capture disease activity, trials may include multiple follow-up times and the method will need to be extended to include latent variables in the mean structure to account for this. In terms of scalability to more complex endpoints, the computational time depends on many things, in particular the number of outcomes, the outcome scale and the number of levels in the ordinal variable. In our case, we find the number of ordinal levels to be the most influential factor in computational time. This is due to the fact that 5 levels in the ordinal variable leads to 10 probability calculations in (3.6), however 3 levels would require the computation of 6 joint probabilities. Consequently, the run time will be substantially increased if there are multiple ordinal levels and decreased if the discrete variables are binary. If the computational time for a particular endpoint is deemed to be too large, then we may reduce the complexity of the endpoint by collapsing the least informative components in to a single binary variable. It must be acknowledged that as we have coded the likelihood, with no package available to do this, the likelihood and probability of response code will have to be tailored specifically to each endpoint. The potential gains in efficiency justify this additional complexity.
We have shown that the latent variable method is a powerful tool in composite endpoint analysis and should be considered as a primary analysis method in a trial using these endpoints. In order for implementation in the general case, where the composite contains any number of continuous and discrete outcomes and to ensure the uptake of the method in clinical trials, we will need to develop a software package. Furthermore, if patients and investigators are to benefit from the efficiency gains, we will need a method to calculate the required sample size in a given trial. We are currently working on addressing these issues to aid in the application of the method. The intuition for the joint probability can be seen below in Figure 10, specifically for the SLE endpoint, where w = 5 and k = 2.
One suggestion in the literature for assessing goodness-of-fit in latent variable models is introduced by Samani and Ganjali (2008) for the case when there is one continuous and one ordinal variable.
This may be extended to allow for two continuous, one ordinal and one binary outcome for application in SLE, as shown below.
As before, let The modified Pearson residuals, taking in to account the correlation between responses are shown below.
A Cholesky decomposition may be used to obtainΣ − 1 2 in (B.1). The covariance between the vector of observed continuous and observed discrete responses is shown below.
y 1 y 2 f Y1,Y2 (y 1 , y 2 )dy 1 dy 2 The Pearson residual is based on the Pearson goodness-of-fit statistics The distribution of the residuals should follow a chi-squared distribution with p degrees of freedom. Comparing the residuals to the chi-squared value allows us to identify observations which the model does not fit well. If there are many observations unexplained by the model then it could indicate a poor choice of model. This may be due to the covariance structureΣ and its assumed distribution. The model may be refitted with various covariance structures and to obtain a model which is found to satisfactorily explain the observed data. If this is not achieved then joint normality of the error terms may be an unreasonable assumption indicating that the latent variable model may not be appropriate. It is possible to fit latent variable models which assume a different multivariate distribution for the error terms, however this is not considered here.

Augmented binary method
The augmented binary model is shown below. The baseline measures for Y i1 and Y i2 are included for comparison, as they are accounted for in the mean structure of the latent variable method.
As one time point is modelled we can use a linear model for Y i1 as shown in (C.1). Note that Y i1 or Y i2 may be chosen as the continuous measure to retain and should always be determined by which is the most informative.
where ε i |T i , y i10 , y i20 ∼ N (0, σ). In this case, the failure time binary indicator will contain information from the remaining three components. F i is set to equal 0 if Y i2 θ 2 , Y i3 is Grade B-E and Y i4 = 0, otherwise the patient is labelled a non-responder in these components and F i1 = 1. F i is modelled using the logistic regression model in (C.2).
Maximum likelihood estimates for the parameters are obtained from fitting models (C.1) and (C.2). The probability of response is shown in (C.3).
As in the latent variable method, (C.3) is used to obtain probability of response estimates for each patient, assuming they were treatedp i1 and not treatedp i1 , which are used to define an odds ratio, risk ratio or risk difference.

Standard binary method
The standard binary method is a logistic regression on the overall responder index, as shown in The odds ratio and standard error estimates can be obtained directly.
Appendix D Relative precision A vs. BV ar(θj ) B V ar(θj ) A θ j : estimated log-odds treatment effect in simulated data j θ: mean log-odds treatment effect over n sim datasetŝ θ low,j ,θ upp,j lower and upper limit of confidence interval for iteration j †V ar[V ar(θ)] =

Scenario
Parameters Investigates Continuous and binary variable driving response Appendix F

Multivariate skew-normal distribution
To test the robustness of the latent variable method to deviations from joint normality of the components, we can generate the data so that the components are drawn from a multivariate skewnormal. The multivariate skew-normal is an extension of the univariate skew-normal distribution introduced by ?. They define it as follows. A random vector Y=(Y 1 , ..., k ) T has k-variate skewnormal distribution, if its density function is where φ k (y; Ω) is the probability density function of the k-variate normal distribution with standardised marginals and correlation matrix Ω. The shape parameter α determines the skewness, where α = 0 reduces the density in (F.1) to the N(0, Ω) density.
Scenarios of interest are shown in Table 7. The first scenario considers when all four components are skewed. Scenarios 2-3 consider different magnitudes of skew in the latent continuous components only. This tests the robustness of the method to the assumption that the observed discrete variables manifest from continuous variables. Scenario 4 is the null case for scenario 3.

Results
The bias, coverage, bias-corrected coverage and power are shown in Table 8 for all four scenarios. In scenarios 1-3, the non-normality introduces bias which results under-coverage. The bias-corrected coverage is close to nominal for all scenarios however the coverage of the latent variable method is nominal in the null case. This is consistent with our findings when the joint normality assumption is satisfied in that bias is introduced in the estimation of the treatment arm, however the magnitude of this bias is much smaller when the assumptions are satisfied. The augmented binary and standard binary methods behave similarly to when the joint normality assumptions are satisfied, which is expected given that the assumptions of those models are violated in both contexts. The latent variable method still offers large power gains over the other methods.  (2017)). Table 10 shows the probability of response in each arm for each of the methods. The findings are consistent with when the assumptions are satisfied. Namely, the latent variable method estimates the probability of response in the control arm well however underestimates the probability of response in the treatment arm. The magnitude of this underestimation is unaffected by the degree of skew or whether the skew is present in the observed continuous components.
The odds ratio treatment effect estimate from each method is shown in Table 11. The latent variable method is biased towards the null, the augmented binary method is biased away from the null. The binary method slightly underestimates the treatment effect in this setting however --all are close to true for the null case.
The median relative precision of the methods are shown in Table 12, with the 10th centile and 90th centile values. These are consistent with our previous findings indicating that the violation of joint normality only affects the bias and not the variance.