{'fraction_non_alphanumeric': 0.02774784482758621, 'fraction_numerical': 0.010876885775862068, 'mean_word_length': 4.866653496641644, 'pattern_counts': {'":': 0, '<': 0, '<?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': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lymphangiogenesis and carcinoma in the uterine cervix: Joint and hierarchical models for random cluster sizes and continuous outcomes

Although the lymphatic system is clearly linked to the metastasis of most human carcinomas, the mechanisms by which lymphangiogenesis occurs in response to the presence of carcinoma remain unclear. Hierarchical models are presented to investigate the properties of lymphatic vessel production in 2997 fields taken from 20 individuals with invasive carcinoma, 21 individuals with cervical intraepithelial neoplasia and 21 controls. Such data demonstrate a high degree of correlation within tumour samples from the same individual. Joint hierarchical models utilising shared random effects are discussed and fitted in a Bayesian framework to allow for the correlation between two key outcome measures: a random cluster size (the number of lymphatic vessels in a tissue sample) and a continuous outcome (vessel size). Results show that invasive carcinoma samples are associated with increased production of smaller and more irregularly-shaped lymphatic vessels and suggest a mechanistic link between carcinoma of the cervix and lymphangiogenesis.


Introduction.
Observational and randomized studies often provide data with a multilevel, or hierarchical, structure, in which repeated data values are available in "clusters" at one level of the hierarchy. Each subunit contributes for data analysis a certain number of observations, which might vary across clusters and which might therefore be regarded as a random variable-a "cluster-specific sample size" or, simply "cluster size". In recent years, consideration has been given to the issue of so-called "informative cluster size", in which the number of observations within a cluster is associated with a study outcome.
This issue is potentially important in many application areas. Prominent amongst them is the field of developmental toxicity, where, for example, a quantitative data obtained from image analysis to describe the number, size and shape of lymphatic vessels in the uterine cervix via measurements of LVD, vessel area and circularity.
Lymphangiogenesis is thought to occur via the sprouting of endothelial cells from existing lymphatic vessels [Alitalo, Tammela and Petrova (2005)]. If this is the case, a subset of smaller lymphatic vessels may be visible in tissue from carcinoma specimens. The structural arrangement of these newly formed vessels in 3D space will influence how functional they are as compared to those found in the normal cervix. This study addresses these issues.
In the context of the more general issue of informative cluster size described above, the methodological challenge arises when further outcome variables-such as the size of each vessel-are associated with the number of vessels at a particular level of observation. For example, in an inverse relationship such as the one described in this paper, large clusters contain vessels that tend to be smaller in magnitude than those in small clusters, yet by their nature the large clusters provide more measurements for analysis. Ignoring the effects of clustering and cluster size may then provide incorrect inferences about the outcome variable, and assessing the extent to which this is true is also one of the objectives of this paper.
Scientific interest in this study therefore lies not only in the number of lymphatic vessels observed, but also in quantitative measures of their appearance. This objective naturally suggests a joint modelling approach. A key requirement is that the model must be flexible enough to allow not only for the correlation among these outcome variables, but also for the fact that they may differ in statistical distribution or data type (e.g., count as opposed to continuous) and at the level of the hierarchy at which they are measured.
In this paper we illustrate the general modelling approach by concentrating on the relationship between LVD and vessel size, adapting for our application methodology that has been successfully applied in fields such as toxicology [Regan and Catalano (1999)]. Our modelling approach relies on the specification of random effects that are common to both outcome variables in the model. These random effects provide a mechanism by which the correlation between the outcomes can be modelled explicitly.
The paper is structured as follows. In Section 2 we describe in detail the data set that motivates this work, and give an exploratory analysis. In Section 3 we fit univariate hierarchical models to each of the outcome variables of interest. In Section 4 we introduce the joint modelling problem and present a bivariate model for lymphatic vessel density and lymphatic vessel area, and Section 5 provides a concluding discussion.  2. Data.
2.1. Study design. The data were collected as part of a study carried out at the Royal Lancaster Infirmary, Lancaster, UK. Tissue biopsies (or "specimens") were taken from 62 individuals. Each specimen was processed into paraffin blocks, which were sectioned at 4 µm, stained and viewed under a microscope, as described by Chapman, Fanshawe and Crick (2013). Within each specimen, areas of interest ("fields") were selected and all lymphatic vessels observed within these fields were used to obtain the outcomes of interest, defined below.
Twenty individuals provided invasive squamous cell carcinoma tissue, while 21 individuals showed premalignant growth classified as cervical intraepithelial neoplasia (CIN), which was additionally subclassified as histological grade 1, 2 or 3 (Table 1). For these two groups of cases, fields were taken from the site of abnormal growth. Additionally, hysterectomy specimens were obtained from 21 controls-defined as individuals with menorrhagia, with no abnormal cervical tissue. For the controls, specimens were available either from one or, more commonly, from two distinct functional regions of the cervix-the ectocervix and the transformation zone.
An average of 8.8 fields were taken per specimen (range 1 to 19 fields per specimen), and across the whole sample an average of 5.5 lymphatic vessels per field provided data for analysis (range 1 to 45 vessels per field). The reasons for apparent differences by group in the number of fields per specimen shown in Table 1 are unrelated to outcome variables.
Two field-level outcomes and two vessel-level outcomes were of primary interest: • Lymphatic vascular density (LVD)-the density of lymphatic vessels visible in a field.
• Percentage lymphatic area (%LA)-the percentage of the total area of a field that is occupied by lymphatic vessels. • Vessel area-the area contained within the lumen of a lymphatic vessel, measured in µm 2 . • Circularity-a measure of the circularity of a lymphatic vessel, lying between 1 (perfectly circular) and 0 (lying in parallel lines across the surface of the field).
Note that as all fields were of the same area, LVD is almost equivalent to the number of lymphatic vessels in a field: a discrepancy would arise only if a convoluted vessel were visible at two or more distinct points on the same field, a scenario that is impossible to detect using the data available and which we consequently ignore. The %LA of a field can be viewed as a combined summary measure of the LVD and the average vessel area of the field, while remaining an important outcome variable in its own right.
A sample size calculation was carried out based on analysis of variance to test for a difference in mean LVD (calculated across all fields, averaging to remove the hierarchical structure) between the three main study groups. Longatto-Filho et al. (2007) and Gao et al. (2006) provide information on LVD in previously conducted studies, although both studies purposively oversampled regions of high LVD. Longatto-Filho et al. (2007) report mean LVD values of 2.6 in the control group (ignoring the distinction between ectocervix and transformation zone), 5.0 for patients with squamous intraepithelial lesions and 17.1 for patients with invasive carcinoma. Because of concerns about the nature of the sampling scheme in these studies, the present study instead assumed more conservative mean LVD values in the three groups of 2.6, 5.0 and 10.0, respectively, which correspond to a between-group standard deviation of 3.8. Based on further results provided by Longatto-Filho et al. (2007), the common within-group standard deviation was assumed to be 7.5, which yields a "difference parameter" [Day and Graham (1989)] of 3.8/7.5 ≈ 0.5. Under these assumptions, a sample size of 25 patients per group has approximately 90% power to detect an overall difference in mean LVD between groups at the 5% level of significance. Incorporating repeated measurements from different fields into the analysis is likely to increase the power substantially.
2.2. Exploratory analysis. Table 2 shows the means and standard deviations of the outcome measures by group. The distributions of three of these variables are illustrated by Figure 1, in which the two control tissue groups are combined, as are the three CIN groups. Vessel area has an extremely positively skewed distribution and is therefore presented on the logarithmic scale.
Two key further considerations guide our approach to analysing this dataset. Firstly, outcomes vary according to the level of the hierarchy at which  they are observed-fields within specimens, and vessels within fields. Figure 2 provides a typical example, and shows the log-vessel area of vessels in the 15 fields taken from the same specimen. In a similar vein, Figure 3 shows the variation in log-vessel area of vessels in the first-numbered fields taken from each specimen. The second major consideration is that the outcome variables themselves are correlated, notably the important variables LVD and vessel area. Figure 4 plots log-vessel area and LVD across all vessels, together with the fitted curve resulting from a generalised additive model fit to the data ignoring the hierarchical structure. The sample correlation between the two is −0.36, with stronger correlation apparent in fields with fewer than fifteen vessels.

Hierarchical models.
3.1. Notation. The data set exhibits a clear three-level hierarchical structure, with vessels nested within fields, which are themselves nested within individuals. For specimen i, field j and vessel k, let Y % ij denote the field-level outcome %LA, Y A ijk the vessel-level outcome vessel area, and Y C ijk the vessellevel outcome circularity. Additionally, let N ij denote the field-level LVD, which can be thought of as representing the cluster-specific sample size. The ranges of the subscripts are i = 1, . . . , n = 62, j = 1, . . . , n i and k = 1, . . . , N ij .
Note that n (the number of individuals) and the n i (the number of fields for individual i) are regarded as fixed, and determined by the study design, whereas N ij is a random variable, an observation that will be explored further in Section 4. Note also that the j subscript enumerates fields in both the ectocervix and transformation zone for the fifteen controls who contribute both of these tissue types. Let x ij denote the tissue type of specimen i and field j, with x ij ≡ x i for CIN and invasive carcinoma specimens. In subsequent regression models results are presented relative to the reference category of control ectocervix. In most cases the three CIN categories were combined for the purposes of model fitting because of the small sample size in the CIN3 group and because differences between the other two categories were small.
The models in this section are conditional on the observed value of N ij , an extremely common approach in the analysis of hierarchical or multilevel data, albeit one that is often made only implicitly. For example, in educational research and in many cluster randomised trials [e.g., Carter (2010)], random cluster sizes are widely discussed, and considerations such as school size and hospital or ward size in a health setting might reasonably be expected to demonstrate an association with outcome measures. In the present study, N ij is prespecified as a key outcome variable to which a priori hypotheses attain.
Different models were formulated for each of the four outcome variables, although each is a variant on the three-level hierarchical model with independent random effects governing within-level correlation [Pinheiro and Bates (2000)].
3.2. Models. For clarity of notation in this section, subscripts are not used to distinguish between corresponding parameters (such as the intercept parameter) for different outcomes. Percentage lymphatic area was modelled untransformed as where a i ∼ N(0, τ 2 ) independently is a specimen-level random effect and Z ij ∼ N(0, σ 2 ) is an error term, where the Z ij are independent of each other and of the a i , that is, fields are assumed to be conditionally independent, given the specimen. Unless otherwise stated, similar independence assump- tions are made on random effect and error terms in the other models described in this section.
The count variable LVD was modelled using a generalised mixed model. An offset of unity was included to prevent zero counts, as any fields containing no vessels would not have been included in the data set, The above formulation can be viewed as part of the general class of models for overdispersion discussed by Dean (1992). In order to test whether this adequately captured the overdispersion in N ij , results were compared to those from a negative binomial model, following Lindsey (1999).
The two vessel-level variables were modelled using an extension of the hierarchy to allow for field-specific random effects. Vessel area was logtransformed and modelled as with b ij ∼ N(0, ν 2 ) independently of the a i and the Z ijk .
The circularity variable Y C ijk required a somewhat different approach owing to the substantive hypothesis relating to this variable. While hypotheses for the other three outcomes all referred to mean differences between tissue types, the question to be answered using the circularity data is whether there is a differential in the structure of lymph cells drawn from specimens belonging to different groups. In particular, it was hypothesised that lymphatic vessels in control tissue would retain a more regular structure than those in case tissue; in this case, control lymphatic vessels might tend to have the appearance of lying parallel to one another. This question is not readily answered using mean circularity, which would instead measure the extent to which lymphatic vessels tend to lie parallel to the plane at which the biopsy is taken, cutting the three-dimensional tissue at an arbitrary angle [Wicksell (1925)].
Instead, it is required to test whether the within-specimen, within-field variation in lymphatic vessel circularity is greater in cases than in controls, which would indicate a less regularly aligned lymphatic vessel network. Thus, using a logit transform to transform the domain of Y C ijk from (0, 1) to

MODELS FOR RANDOM CLUSTER SIZES AND CONTINUOUS OUTCOMES 11
(−∞, ∞), the model is of the form where random effects are specified as above, except that b ij ∼ N(0, δ x ij ν 2 ), with δ 4 (corresponding to invasive carcinoma) set to unity for identifiability, and we require to test the hypothesis that all of the δ-parameters are equal against a general alternative.
Models were fitted using functions in the nlme and MASS packages in R [Pinheiro et al. (2008), Venables and Ripley (2002), R Development Core Team (2008)].

3.3.
Results. The results of fitting models (1)-(4) are summarised in Table 3. For all outcomes, there was statistically significant clustering within specimens and also within fields (nonzero random effect variance parameters). Intra-cluster correlation coefficients can be estimated from the variance parameter estimates in Table 3; for example, the within-specimen clustering effect for %LA is estimated as 1.20/(1.20 + 8.63) = 0.12. For the two field-level outcomes, clustering effects are higher at the field-specific level of the hierarchy than within specimens. Standard errors of most fixed effect estimates are increased by a factor of around two compared to corresponding models that make no allowance for the hierarchical data structure, and point estimates remain similar. No improvement on the Poisson random effects model was seen using the negative binomial distribution (in terms of change in the log-likelihood), and so the Poisson formulation was retained.
Compared to control ectocervix tissue, there is clear evidence that all other tissue types have greater LVD and smaller vessel area. The difference is especially marked for invasive carcinoma tissue, for which there is an estimated 3.7-fold increase in LVD and 3.8-fold reduction in average vessel area compared to control ectocervix. As LVD and vessel area show opposing trends, as a result of the correlation between them, there are no differences between normal ectocervix, CIN and invasive carcinoma tissue for %LA, although there is evidence that %LA is higher in normal transformation zone tissue.
There is evidence that vessels tend to be more circular in control transformation zone and CIN than in the other two tissue types, but also evidence that the within-field variance of circularity measurements is higher in invasive carcinoma specimens than in control ectocervix and CIN (as estimates of the relevant δ-parameters are significantly less than unity). Moreover, analysis of variance comparing the fit of this model with the special case in which all δ-parameters are constrained to be equal suggests significantly improved fit of the more general model (p = 0.001). Allowing a separate δ-parameter for each of the three CIN categories, however, did not greatly improve the fit (p = 0.07), so the more parsimonious model is reported here. It should be noted that the study was not designed to investigate differences between the three CIN grades. The final model thus suggests greater variation in circularity amongst fields taken from invasive carcinoma tissue than for control ectocervix and CIN.

Introduction.
The key issue remaining to be addressed relates to the assumption made implicitly in hierarchical modelling such as in Section 3 that the cluster-specific sample sizes are fixed. In the present study it is reasonable to regard the total sample size, n, and the number of fields per individual, n i , to be fixed by design, but the number of vessels per field, n ij , would be more accurately regarded as the realisation of a random variable N ij , which is the definition of the LVD outcome variable.
Fitting models that condition on the observed value of this random variable may lead to incorrect inferences about the relationship between study group and other outcomes of interest. For example, there may be underlying tissue-specific characteristics that are associated with both an increase in LVD and a reduction in vessel size; indeed, it is biologically plausible that this should be the case in invasive carcinoma tissue. These effects may be masked by an analysis that is based on conditioning on LVD, and the conditional models may give inappropriate inferences about the relationship between tissue type and vessel size. 4.2. Models. One possible analysis strategy factorises the joint likelihood of Y and N into conditional and marginal components, that is, [Y, N |X] = [Y |N, X][N |X], where the notation "[·]" means "distribution of". This has the appealing property of allowing a simple marginal analysis of N |X, but raises the question of whether the models for Y |N, X considered in Section 3 are satisfactory even as conditional models, as they do not explicitly model the correlation between N and Y .
As a possible solution, Catalano and Ryan (1992) discuss bivariate models in which a function of N enters the expression for [Y |N, X] directly as a covariate. This approach is also adopted by Panageas et al. (2007) in the context of outcome of surgery when the cluster size is the number of patients treated for a given surgeon.
Geometrical considerations and Figure 4 suggest that an appropriate choice in the present study might be N −1 for the following reasons. In 2D, lymphatic vessels are approximately circular in appearance. The maximum total area of m 2 circles of equal radius packed into a square of area a 2 is πa 2 /4m 2 , and so in two dimensions the proportion of the field area that is filled by vessels might be expected to be proportional to the reciprocal N . However, for reasons discussed above, there are drawbacks to a conditional model, which in any case answers a different research question to the one of primary interest, and this simple approach also fails in itself to take account of the stochasticity of N .
An alternative is to model the joint distribution of Y A ijk and N ij directly. The correlation can be modelled using a linked random effect approach: where b ij and Z ijk are as previously and, independently of b ij and Z ijk , This formulation closely follows the general three-outcome model for one continuous, one binary and one count variable discussed by Catalano and  (1992) and Gueorguieva and Agresti (2001), and developed further by Gueorguieva (2005). The shared random effect structure considered by these authors is very similar to (5), although in place of the Poisson model for N ij they instead use a continuation ratio probit model, which in our notation would take the form

Ryan
introducing a potentially large number of additional parameters δ h . The work of Gueorguieva (2005) extends the model of Dunson, Chen and Harry (2003), who provide an integral expression for the correlation between the two outcomes conditional on the random effects. An advantage of this joint approach over the conditional one is that it gives an estimate of the direct effect of tissue type on each outcome variable, similar to that used in joint longitudinal and survival modelling [Ibrahim, Chu and Chen (2010)].
It is convenient to use a Bayesian framework to fit model (5). Priors were specified as follows: for each α and β parameter, N(0, 10 −6 ); for each precision parameter, Γ(10 −3 , 10 −3 ); λ ∼ Unif(−10, 10); ρ ∼ Unif(−0.95, 0.95). The prior for ρ was chosen with two considerations in mind: to allow that the outcomes might be strongly, but not perfectly, correlated, and to ensure that draws from the distribution of ρ do not allow Σ to become singular. The model was fitted in WinBUGS v1.4 [Lunn et al. (2000)], with a burn-in period of 50,000 iterations. Posterior estimates were obtained from a further 50,000 iterations, with a thinning factor of 20. In additional analyses to check the sensitivity to the choice of priors, the prior distributional forms and/or numerical values of the hyperparameters were varied (although always remaining "vague", in the sense of having high variance).
For comparison, this model was also fitted with ρ set to zero, equivalent to fitting univariate models to each of the outcomes in the Bayesian framework. Table 4 shows parameter estimates (median of posterior distribution and 95% credible interval) from model (5). Fitting the model with ρ set to zero gives near-identical results to Table 3, the only difference being the LVD Transformation zone parameter, for which the estimate in the Bayesian model was 2.37 [1.96, 2.84]. This allows direct comparison of the univariate and joint results. After thinning, the autocorrelation of posterior samples for all parameters was negligible, and there was no material difference in the results according to the choice of priors.

Results.
Point estimates for the β parameters are similar to those in Table 3, but in many cases credible intervals are somewhat wider. The negative estimate of λ N was expected and is due to the negative correlation between LVD and vessel area. Similarly, the estimate of ρ is negative, and suggests high correlation between the individual-level random effects. The posterior distributions of most parameters, including the components of the random effects b, are approximately symmetric. An exception is ρ, by virtue of being bounded through the support of the prior by ±0.95, which therefore has a positively skewed posterior distribution. Figure 5 shows the relationship between the fitted field-level random effects from model (3) and LVD. This provides a possible reason why the point estimates of the joint model are largely unchanged compared to those from the single-model, as the relationship is strikingly similar to that between log-vessel area and LVD (Figure 4). In the univariate model the correlation between the cluster size and the outcome is implicitly accounted for through the random effects distribution, even though no such association has been specified in the formulation of the model.

5.
Discussion. The analysis presented provides clear evidence of differences in properties of lymphatic vessels according to presence of CIN or invasive carcinoma in the cervix. Relative to the other tissue types considered, tissue taken from individuals with invasive carcinoma have lymphatic vessels that are greater in number, smaller in size and less regular in shape. In addition, this study found a difference between the ectocervix and transformation zone of the control cervix, whereas no difference was found in LVD between control transformation zone and CIN groups. The main scientific conclusions were relatively unaffected by the decision to model the data in a joint as opposed to a univariate framework.
These findings are consistent with a model that asserts that lymphangiogenesis occurs when the cervix undergoes eversion at puberty, to which all experimental groups would be subject, but also suggest that there is a separate lymphangiogenenic episode in the squamous carcinoma group. Further investigation is needed to ascertain whether the expression of genes encoding growth factors that may cause this additional lymphangiogenesis occurs late in CIN or after the progression to squamous cell carcinoma. Also, the average area of vessels in the CIN group was unexpectedly smaller than those in the transformation zone, a finding not explained by this model of lymphangiogenesis, confirming the need for further studies. In the carcinoma group, vessels were considerably smaller than in all other groups, supporting the hypothesis of newer formation, and there was some evidence that they are morphologically and therefore functionally different.
The approach to statistical analysis was guided by the multilevel structure of the data set, in which it was required to account for high within-specimen correlation in the outcome variables measured. The negative correlation between LVD and vessel area suggested a joint model for these two outcomes would provide a measure of the direct effect of tissue type on the two outcomes. As one of these variables was equal to the cluster-level sample size, a different analysis strategy was required than that typically used, for example, for two continuous outcome variables.
Our method used shared random effects in a bivariate Normal-Poisson framework to allow for this correlation. We used a Bayesian framework, al-though similar models have been considered utilising direct likelihood maximisation [Gueorguieva (2005)]. Allowing for the randomness in cluster size tended to increase the standard errors of parameter estimates, but there were no substantial changes in the point estimates themselves.
This finding agrees with the results of Neuhaus and McCulloch (2011), who demonstrate that for linear mixed models with random intercept only, estimators of covariate effects are consistent, and are estimated equally efficiently, even when cluster size is ignored. They point out that the exception is the fixed intercept parameter (α), which in any case is rarely of substantive interest. However, they assume that cluster sizes do not depend on covariates, which is not plausible in our application in that cluster size itself is an outcome measure, with strong evidence of an association with tissue type.
Additionally, the data considered here did not warrant using random slope models, which have been used successfully elsewhere when continuous covariates were available [e.g., Dunson, Chen and Harry (2003)]. In general, the area of informative cluster size has received little attention in the statistics literature, and has only been applied to a limited number of application areas, which is perhaps surprising given the rapidly increasing level of research in the area of cluster randomised trials.
The joint modelling framework considered could in principal be extended to the multivariate case. In our analysis this was not necessary, as there was no reason why circularity should be associated with the number or size of vessels. To achieve this, similar models might be considered that use several shared random effect terms to induce correlation between the random variables. Catalano and Ryan (1992) provide further details.
This study did not attempt to distinguish between the cancer grade or stage of samples analysed, and also makes no link between properties of lymphatic vessels and metastasis or prognosis. These are both limitations and possible future research directions, and it would be a simple extension of the joint framework set out here to analyse data of this type. In addition, in this study it was not possible to measure either the spatial distribution of lymphatic vessels within a specimen or changes in lymph structure relative to distance from the tumour site, both of which may provide extra insight.