Multilevel Latent Class Modelling of Colorectal Cancer Survival Status at Three Years and Socioeconomic Background Whilst Incorporating Stage of Disease

Previous studies investigating survival from colorectal cancer have typically considered potential confounders to include stage of disease. Stage however may lie on the causal path and statistical adjustment with stage as a confounder can then introduce bias known as the reversal paradox. Classification of stage may also be imprecise and incomplete. Modelling using Latent Class Analysis (LCA) may minimise bias by including covariates on the causal path as ‘class predictors’ and by accommodating uncertainty associated with confounder values explicitly via the latent class part of the model. We construct multilevel latent class models to allow for the multilevel structure of the data: patients nested within NHS Trusts. We use a dataset of patients in a large UK regional population diagnosed with colorectal cancer between 1998 and 2004. Death within three years is the outcome. The optimum number of latent classes at patient and Trust level is determined with reference to likelihood-based modelfit criteria. The three-patient five-Trust class multilevel LCA model was chosen. Patient classes were identified as good, reasonable or poor prognosis groups. The impact of stage differed across the patient classes. Socioeconomic background and older age were clearly associated with increased odds of death in all patient classes. Females had significantly decreased odds of death compared with males in the good prognosis class. The five Trust classes identified outlying Trusts, indicating that the standard multilevel model would not have been sufficient to model these data.


Background
Many factors influence survival from colorectal cancer, including diagnosis or treatment centre (see Kee et al., (1999) on the influence of hospital and clinician workload; McArdle & Hole (2004) on volume and specialisation; and Borowski et al. (2010) for a volume-outcome analysis), stage at diagnosis (see Woodman, Gibbs, Scott, Haboubi, & Collins (2001) on differences in stage at presentation; and Ciccolallo et al. (2005) on the role of stage and surgery), and associated risk factors such as socioeconomic background, age at diagnosis and sex (see Morris et al. (2011) for thirty-day postoperative mortality; Downing et al. (2013) for early mortality; Smith et al. (2006) on the impact of social deprivation; Widdison, Barnett, & Betambeau (2011) on age; and Hendifar et al. (2009) on gender disparities).Findings exploring the potential impact of socioeconomic background (SEB) vary, with some studies showing an impact of poor SEB on decreased colorectal cancer survival (Downing et al., 2013;Morris et al., 2011), whilst others do not find this association (Nur et al., 2008;Smith et al., 2006).We investigate the relationship between survival status from colorectal cancer and SEB, while accounting for other factors that may affect this relationship using a novel statistical approach that can account for both genuine confounding and what is effectively mediation, often mistakenly treated as confounding but which may in fact bias the estimated impact of SEB when adjusted for in standard regression models.
For instance, previous studies may have included stage of disease at diagnosis as a potential confounder to the relationship between known potential risk factors and survival from colorectal cancer.However, a higher level of deprivation may provoke late presentation, perhaps due to the refusal of a screening invitation (Whynes, Frew, Manghan, Scholefield, & Hardcastle, 2003), which may result in a more advanced stage at diagnosis (Ionescu, Carey, Tait, & Steele, 1998).SEB therefore causally precedes stage at diagnosis and consequently stage does not qualify as a genuine confounder; it is a mediating factor.Bias may be introduced by the statistical adjustment for mediators on the causal path (Kirkwood & Sterne, 2003), termed the reversal paradox (Stigler, 1999), which may be a serious problem in epidemiology (Hernández-Díaz, Schisterman, & Hernán, 2006;Tu, West, Ellison, & Gilthorpe, 2005).Figure 1 shows a theorised diagram of causality using a Directed Acyclic Graph (DAG) (Pearl, 2000) to determine which variables in our study are confounders, proxies for confounding, competing exposures, or mediating variables that lie on the causal pathway., 2008).Classification may also be imprecise as patients may be classified incorrectly due to the variable quality of pathology (Quirke & Morris, 2007) or be "understaged" (i.e.incorrectly assigned a earlier stage at diagnosis due to unidentified lymph node metastases) (Morris, Maughan, Forman, & Quirke, 2007).Statistical analyses using regression modelling may yield biased results where model covariates have measurement error (Greenwood, 2012) or missing values (Carroll, Ruppert, Stefanski, & Crainiceanu, 2006;Fuller, 1987), and this bias is exacerbated when considering product interaction terms (Greenwood, Gilthorpe, & Cade, 2006).Models that incorporate staging data may therefore introduce bias due to the variable quality and completeness of pathology.

Modelling Approaches
Regression modelling (Normand et al., 2005) is often extended to a multilevel framework in order to incorporate differences across diagnosis or treatment centres (Leyland & Goldstein, 2001;Leyland & Groenewegen, 2003), such as NHS Trusts.This approach however assumes that a study sample is homogeneous at every level, the same model would be applied to all members of the sample, and the effects of covariates would thus be the same throughout.In a multilevel model, variation amongst the intercepts and slopes is assumed to be normally distributed and independent of the variation in the individual measurements.These assumptions may not be valid in observational health data; patients or Trusts are unlikely to be homogeneous where patients have not been randomly selected for inclusion and Trusts have not been randomly distributed geographically.As such, it may be inappropriate to apply one model to all individuals, and the effect of covariates therefore cannot be assumed to be the same throughout the sample.
Latent Class Analysis (LCA) could be used, and this also extends to a multilevel framework.Downing, Harrison, West, Forman and Gilthorpe (2010) incorporated stage at diagnosis as a "class predictor" in the LCA rather than as a covariate and found improvements in model fit using multilevel LCA (MLLCA) in comparison to singlelevel logistic regression modelling when studying risk factors related to breast cancer survival status.MLLCA is therefore proposed here to illustrate an original application in an area where its utility may be overlooked.It is important to consider alternative approaches to match the context of the data; we advocate an improved approach to analysis in cancer research data.We construct multilevel latent class models to identify subtypes of patients and NHS Trusts, simultaneously, to model how patients may vary and how NHS Trusts may differ, based on survival status.We compare the MLLCA approach with standard multilevel models (MLMs) to examine improvements in model fit and model interpretation and hope to demonstrate the utility of the latent class approach.

The Colorectal Cancer Dataset
The Northern and Yorkshire Cancer Registry and Information Service (NYCRIS) database was used to identify cases of colorectal cancer (ICD-10 (World Health Organisation, 2005) codes C18, C19 and C20) diagnosed between 1998 and 2004, where the patient was resident in the Northern and Yorkshire regions.A description of the data extraction and exclusions is available in a previously published study (Gilthorpe, Harrison, Downing, Forman, & West, 2011).The outcome was whether or not the patient survived at three years following diagnosis, as this is clinically meaningful and facilitates ready comparison with other studies.Following exclusions 24640 records were available for analysis.Data include information on age at diagnosis, sex and SEB (using the Townsend Index (Townsend, Beattie, & Phillimore, 1987) recorded at the 2001 census, stage at diagnosis (using the Dukes classification (Dukes, 1949)), the ICD-10 diagnosis code for the tumour, its laterality (position in the body), and whether or not the patient was treated curatively.The diagnostic centre was defined as the NHS Trust where the latest staging took place; 19 Trusts were identified in the NYCRIS geographical area.
The colorectal cancer data are hierarchical since different groups of patients attend different diagnostic centres, dependent on factors such as their area of residence; patients are clustered within NHS Trusts and there will be variation at both patient and Trust level.A multilevel modelling framework would therefore seem appropriate.It is also important to account for stage at diagnosis, but this cannot be included as a covariate alongside SEB without the risk of introducing bias due to the reversal paradox.

Latent Class Analysis
LCA is also known as discrete latent variable modelling, or mixture modelling (Goodman, 1974;Magidson & Vermunt, 2004).In LCA, a number of latent classes, or subgroups, are identified, the optimum choice of which is selected by the researcher (usually informed by log-likelihood statistics).Units of analysis are assigned to a latent class based on similarities in their characteristics and latent classes are therefore homogeneous, with similar effects of each covariate on units in the same latent class, though covariate effects may differ across the classes.The relationship between outcome and associated risk factors can thus be determined within each latent class, rather than over all observations.
LCA has the utility to model covariates as "class predictors".This may be either in addition to or instead of their inclusion as standard covariates along with the main exposure under investigation.For confounders that are also potential effect modifiers (i.e. they exhibit an interaction with the main exposure), modelling these variables as class predictors yields an implicit interaction, since the outcome-exposure relationship may vary across latent classes.This averts the need to include an explicit confounder-exposure product term in the standard part of the regression model, which would otherwise exacerbate any bias introduced if the confounder is measured with error or has missing values.Modelling effect modification this way minimises bias; uncertainty associated with confounder values is explicitly accommodated via the latent class part of the model.

Confounding and Mediation
If an alleged confounder lies on the causal path between exposure and outcome, it is a mediator, and its statistical adjustment in the standard regression model introduces bias; it is then wise to discard the mediator as a model covariate.This does not preclude the mediator becoming a "class predictor", though one then has to ensure there is no remaining implicit bias.Modelling a mediator as class predictor yields the potential for implicit interaction, as before, where the exposure-outcome relationship may vary across latent classes.The exposure may thus cause the mediator, which in turn part determines the latent class structure, within which the exposure-outcome relationship may vary.Circularity thus arises in the causal interplay of exposure, mediator and outcome.This can be avoided if the outcome-exposure relationship is not allowed to vary across latent classes.In such instances, only the intercept varies across each latent class, not the exposure-outcome slope.Although the causal circularity is avoided, this may not avoid some degree of residual bias due to the reversal paradox, as the exposure-outcome relationship is unlikely to be independent of within-class intercepts, which effectively are 'adjusted' by the consideration of the mediator as a class predictor.We nevertheless explore the notion that variables which lie on the causal path between exposure and outcome may be considered as class predictors instead of being incorrectly adjusted for as alleged confounders within the standard regression model.

Multilevel Latent Class Analysis
MLLCA is an extension of LCA with latent classes determined at more than one level, where classes at the lower level are based on similarities in characteristics (Skrondal & Rabe-Hesketh, 2005).Latent classes at the lower level are homogeneous, while those at the upper level can be homogeneous or heterogeneous, dependent on model specification and research question.An optimum model is sought for all classes at all levels simultaneously.
Covariates can be included at any level and, as with single-level LCA, their effect is the same within each latent class but may differ across the classes (if deemed appropriate).If intercepts and slopes are fixed within classes at all upper levels, no distributional assumptions are required.As within single-level LCA, covariates can be modelled separately to the main association under investigation, as 'class predictors'.
We use MLLCA to analyse the colorectal cancer data.Ultimately, no assumption of normality is made at the upper level, though initially a continuous latent variable for the upper level is adopted (as an approximation) while the latent class structure is explored for the lower level.Once the lower-level optimum number of classes is determined, the upper level latent variable is switched to categorical and the optimum upper-level latent class structure is determined.Stage at diagnosis is considered as a class predictor with the SEB-survival relationship held constant across the patient-level latent classes, thereby minimising the effect of the reversal paradox and potential bias introduced due to measurement error in the stage variable.Trust classes are homogeneous with respect to both patient outcome and its relationship with model covariates.This places the focus on patients and allows us to determine what kind of patient is potentially susceptible to the adverse impact of SEB in terms of their cancer survival status (i.e.we determine patient casemix characteristics in relation to outcome).An alternative approach of grouping Trusts according to differences in characteristics is discussed by Gilthorpe et al. (2011), where differences in survival at the Trust level may be as a result of underlying differences in Trust performance, rather than patient casemix.

The Modelling
We use both likelihood-based model-fit criteria and a graphical method to determine the optimum number of latent classes at each level; we consider both the Bayesian Information Criterion (BIC) (Schwarz, 1978) (for reasons of parsimony) and the change in log likelihood (LL).We also examine and report classification error (CE), but do not use it to inform our model choice.CE reflects the proportion of misclassified observations (at each level separately) when comparing the modal and probabilistic assignment to classes; a lower CE signifies that the latent classes are more "real", i.e. observations are almost entirely assigned to single classes.Models include adjustment for age at diagnosis, sex and SEB.An age-squared term is also included as age was found to have a non-linear relationship with survival; the inclusion of age-squared allows for an adjustment to the linear effect of age.Stage at diagnosis is a 'class predictor', and the ICD-10 diagnosis code, laterality and whether treatment is curative or not are 'inactive' covariates at the patient level, i.e. not used to estimate the model, though used to partition the findings by these variables for descriptive purposes.We generate 200 bootstrapped datasets and analyse each similarly in order to generate 95% confidence intervals (CIs).The software Stata (StataCorp, 2011) was used for data manipulation, summary statistics, tabulation and charts, while the statistical software LatentGOLD (Vermunt & Magidson, 2005) was used for the latent class analyses.Table 1 shows the results of the MLM analysis.Overall, 12 708 patients (51.6%) died within three years.The reference group comprised males aged 71.5 years (the mean age), diagnosed with Stage A colorectal cancer and with a zero Townsend score.Substantial and statistically significant associations were found between increasing deprivation and increased odds of death (OR = 1.18, 95% CI = 1.15 to 1.21 per SD increase in Townsend score); between female gender and decreased odds of death (OR = 0.87, 95% CI = 0.83 to 0.92); between increasing age and increased odds of death (OR = 1.31, 95% CI = 1.30 to 1.33 per 5-year increase in age); and between increasing age squared and increased odds of death (OR = 1.006, 95% CI = 1.005 to 1.007).

MLM Analysis and LCA Approach
With a continuous latent variable at the upper level, the MLLCA approach suggests that three patient classes are optimum by both the BIC and change in LL.Table 2 summarises the model-fit criteria for the multilevel latent-class models on switching the upper level latent variable to a categorical to determine the optimum upper-level latent class structure, and shows the optimum models identified by each criterion.Table 2 shows that one Trust class is optimum by the BIC.More than one class at the Trust level is required to explain Trust differences however, therefore we consider also the -2LL plot shown in Figure 2 , 2004).There is little indication that either the type or position of the tumour is associated with survival status, as the proportions are broadly similar across the patient classes.
Across all patient classes, SEB was clearly associated with increased odds of death (OR = 1.33, 95% CI = 1.26 to 1.41).In the good prognosis patient class, females had significantly decreased odds of death compared with males (OR = 0.59, 95% CI = 0.40 to 0.87); in the reasonable and poor prognosis classes the association was less clear Interval; mean Townsend score over all classes -0.04; mean age over all classes 71.5 years.CIs from bootstrapping calculated using percentiles.
Table 4 summarises the chosen model for the Trust classes, where Trusts were apportioned into one of five groups.According to modal assignment, class 1 contained seven Trusts (37.3% of patients); class 2 contained five Trusts (26.9% of patients); class 3 contained three Trusts (14.3% of patients); and classes 4 and 5 contained two Trusts each (11.1% and 10.4% of patients respectively).According to probabilistic assignment, the prevalence rates ranged from 49.6% of patients dying within three years (in class 4) to 54.5% (in class 5).
The remainder of the results pertain to probabilistic assignment.SEB differed somewhat across the Trust classes, with the highest value seen in class 3 (mean SEB = 0.38), indicating that Trusts in this class receive patients on average from more deprived areas.In contrast, the lowest values of mean SEB are seen in classes 4 and 5 (mean SEB = -0.39 and -0.45 respectively), indicating that Trusts in these classes receive patients on average from more affluent areas.The mean age of patients remains fairly constant across the classes, ranging from 71.2 years (in class 4) to 71.8 years (in class 2).The proportion of females also remains fairly constant across the classes, ranging from 43.4% (in class 4) to 44.6% (in class 5).No substantial trend is seen across the Trust classes by stage, although class 5 contains slightly more patients with missing values for stage (15.7%)compared with the other classes.The fewest patients received curative treatment in Trust class 5 (81.7% compared with 84.7%, the highest proportion, in class 2), which perhaps indicates that the two Trusts in this class are not treating as many early-stage patients as other Trusts.There are some modest differences seen across the Trust classes in type of tumour: class 5 has the highest proportion of colon tumours (61.8%) and the lowest proportion of rectosigmoid junction tumours (8.9%); while class 3 has the highest proportion of rectum tumours (31.8%).There are also some small differences seen in the laterality of the tumour: class 3 has the highest proportion of tumours on the left side (67.2%); classes 1 and 2 have the highest proportion of tumours on the right side (27.7%); and class 5 has the highest proportion of tumours across both sides (10.8%).

Findings
The MLM analysis found sizeable and significant associations between increasing deprivation and increased odds of death, between being female and decreased odds of death, and between older age and increased odds of death.The MLLCA categorised patients into three latent classes, labelled as good, poor and reasonable prognosis (relating to stage at diagnosis) and in all classes, the impact of the covariates was found to agree with that seen in the MLM analysis.There were some differences across the prognosis groups, with the good prognosis group showing clearly decreased odds of death for females compared with males; and the reasonable prognosis group showing a greater impact of older age.The association between SEB and survival status was clear, with higher deprivation associated with increased odds of death in all patient classes (as the impact of SEB was deliberately held constant across the classes, i.e. there was no stage-SEB interaction, because this would otherwise have introduced circularity in the causal relationships amongst SEB, stage at diagnosis, and the latent classes within which the SEB-survival relationship varied).As discussed, previous findings into the association of SEB with survival from colorectal cancer have been seen to vary and this may in part depend upon whether these studies have undertaken appropriate statistical adjustment for alleged confounders, or introduced bias due to the reversal paradox.
MLLCA has considerable utility to account for issues of structure, non-homogeneity, inferred causality, missing values and measurement error, while improving upon the MLM approach by producing models that are better fitting to the data and provide an enhanced interpretation of the data.It allows for the hierarchical structure of the data without imposing any distributional assumptions.It accounts for non-homogeneity at all levels by categorising both patients and Trusts into latent classes and allowing the relationship between survival status and associated risk factors to vary across these classes (where appropriate), rather than modelling the relationship over all patients and across all Trusts.By modelling stage at diagnosis as a class predictor, bias due to the reversal paradox is minimised.As patient classes depend on stage, we investigate covariate-outcome associations within sub-categories of stage without introducing product interaction terms, thereby minimising bias due to measurement error.We take account of incomplete data within stage by categorising missing values and the modelling assigns patients with missing stage values to the most appropriate patient class according to similarities in their other characteristics compared with other patients.
The determination of Trust classes in the MLLCA was not straightforward.The best fit according to the more parsimonious measure of the BIC was one Trust class but this would not be sufficient to describe fully the natural structure of the data.Model fit according to the LL improved as the number of Trust classes increased, surpassing the fit of the standard multilevel model when considering only two Trust classes; the -2LL plot showed that the model fit continued to improve up to around five Trust classes (Figure 2).These five classes identified outlying Trusts, with predominantly three Trusts in Trust class 3 and two Trusts in each of Trust classes 4 and 5, indicating that the normal approximation at the upper level would not be ideal.This confirms our suspicion that the standard multilevel model would not be the most appropriate to model these data.
The five Trust classes differ only due to patient survival status and the relationship between survival and the covariates, potentially enabling us to highlight differences in patient care (e.g.treatment pathways or hospital characteristics, such as size or speciality) that might explain the differences and so be worthy of further investigation.It should be noted that the differences seen across the Trust classes are not statistically significant.Nevertheless, their inclusion was necessary in order to account for variability at the Trust level and so to model best the corresponding patient classes.No Trust level covariates were available for inclusion in the modelling, meaning that this potential aspect of investigation could not be addressed in this study.

Limitations
Although we minimise any potential bias that arises due to the reversal paradox by including stage at diagnosis as a class predictor, and by holding constant the SEB-survival relationship across classes, this bias may not be entirely eradicated and it is unclear how much bias could remain.No methodology can completely eliminate all bias and the use of MLLCA will have reduced the risk of bias in comparison with standard MLM techniques.
Although interest may lie in establishing treatment centre characteristics that may have an impact on survival from colorectal cancer, we have modelled diagnostic centre at the upper level.This allowed us to include all patients regardless of whether or not they received treatment.We modelled Trust of diagnosis at this level to minimise error that could be introduced by patients receiving treatment at different hospitals during their care, as a higher proportion of patients receiving treatment were treated within the same Trust at diagnosis (90.2%), than within the same hospital at diagnosis (81.7%).
We have included SEB at the patient level to simplify the analysis, though it is derived at the small area level.Individual measures of deprivation are rarely available, especially when using routine data.Extrapolating area-based findings to individuals, however, can lead to the ecological fallacy (Robinson, 1950).An additional upper level could be introduced relating to patients' area of residence and this would be cross-classified with Trusts, since patients attending hospitals in one Trust may be resident of different areas.We could also consider survival as a continuous measure, as within Cox proportional hazard modelling.Both these extensions could be accommodated within a MLLCA framework, though alternative software such as WinBUGS (Lunn, Thomas, Best, & Spiegelhalter, 2000) or MPlus (L.K. Muthén & B. O. Muthén, 1998-2011) would then be required (although cross-classified extensions are still under development for MPlus).For simplicity of illustration of the methodology, however, we did not pursue these options within this study.
For this paper, we chose to focus on the utility and interpretability of an alternative modelling approach that has the potential to model appropriately both confounders and mediators while preserving the clinical message.Simulation studies may be beneficial to gain insight into the sensitivity of the data to model choice, and these could be considered as an extension to this study.

Conclusion
The MLLCA modelling approach illustrated a better fit to the data and showed new insights that were not previously apparent using the MLM approach.The impact of covariates on survival status differed across latent classes defined by stage at diagnosis.By tailoring treatments and pathways according to patients' profiles, there might be opportunities in the future to optimise patient care.This analytical strategy has prognostic utility to inform health service providers of disparities within patient care.

Figure 1 .
Figure 1.Directed Acyclic Graph (DAG) showing the inferred causal relationships amongst all available variables at the population level , which shows model fit improving as the number of Trust classes increases.The standard MLM showed a LL of -11 985 which is surpassed by using two Trust classes.In order to model fully Trust variability and to improve patient class estimates we choose the model with five Trust classes, which lies at the point where there is little further improvement in model fit.

Figure 2 .
Figure 2. -2LL plot used to determine the optimum number of Trust classes in the MLLCA approach

Table 1 .
Results from MLM analysis (multilevel logistic regression): odds of death within 3 years Odds Ratio, CI-Confidence Interval, SD-Standard Deviation; the reference group consisted of males of mean age, diagnosed with Stage A colorectal cancer and with a zero Townsend score; LL = -11 985.

Table 2 .
Model-fit criteria for the three-patient multilevel latent-class models with a categorical upper level latent LL-Log likelihood; BIC-Bayesian Information Criterion; CE-Classification Error.

Table 3 .
Results for the patient classes in the three-patient, five-Trust-class multilevel latent-class model Confidence Interval, SD-Standard Deviation; the reference group comprised males, aged 71.5 years, classified as Stage A at diagnosis and attributed a Townsend score of zero; the Wald p-value indicates levels of statistical significance for differences in effect across the classes.CIs from bootstrapping calculated using percentiles.Table3summarises the patient classes from the chosen three-patient five-Trust-class MLLCA model, where patients were apportioned into one of three groups, labeled post-hoc as: good prognosis, reasonable prognosis, or poor prognosis.The good prognosis class contained 38.2% of cases of which 9.4% died within three years, compared with the reasonable prognosis class with 27.6% of cases of which 58.3% died within three years, and the poor prognosis class with 34.2% of cases of which 93.2% died within three years.The profile of stage differed across the patient classes.The good prognosis class corresponds to early-stage diagnosis with 70.8% of the stage A/B patients.The reasonable prognosis class corresponds to mid-stage diagnosis with 59.1% of the stage B/C patients, and a large proportion of patients with missing values for stage (30.5%).The poor prognosis class corresponds to late-stage diagnosis with 82.6% of the stage C/D patients.The good prognosis class contains the highest proportion of patients treated curatively (98.8%), which may be partly due to their stage at diagnosis as early-stage patients commonly receive curative, instead of palliative, treatment (National Institute for Clinical Excellence