The New England Journal of Statistics in Data Science

Networks are collections of interconnected entities, e.g. social networks of communicating actors, ecological networks of flora and fauna commensalism, and computer networks.[31, 26] Network-graphs, or simply graphs, are mathematical objects consisting of a vertex set, e.g., one vertex per network entity, and an edge set, e.g., a set of pairs of vertexes involved in a network connection, that represent the arrangements of pairwise relationships in the network.[9]

Methodology developed in the fields of social networks, network science, and graph theory provides for the analysis of relational data generated from a variety of measurements across scientific disciplines.[22, 8] Community detection is the process of identifying exceptionally dense subnetworks of mutually well-connected network entities, known as communities, that often have functional meaning in the network.[12] Notable approaches to community detection include the clique percolation method [24], spectral partitioning [5], degree-corrected stochastic block models [18], modularity optimization [21], and multi-slice network community detection [20]. These approaches are designed for the unsupervised partitioning of the vertex set of a graph into unusually cohesive subsets of vertexes, and with varied applications in sociology [33], computer architecture [13], and biology [2], illustrate that myriad methodology in community detection procedures are applied broadly in scientific research.

Community detection procedures commonly integrate network connectivity exclusively, without regard for other quantities of interest, e.g., auxiliary measurements on vertexes.[10]. We refer the reader to recent expositions on the state of the science of community detection in networks including [16, 30, 34, 27] that detail existing methods with respect to certain applications. It is essential to note that constraints are not commonly imposed on the community structure in networks by existing methods.

We are motivated to partition a nation-wide network of hospitals into subnetworks of hospitals that (i) exhibit a high level of within-group patient sharing as quantified by the network modularity quality function and (ii) consist of at least one hospital that has hosted a minimum number of implantable cardioverter defibrillator (ICD) surgeries to define as health care communities (HCCs) that provide a comparable level of cardiac care. We utilize data acquired from health insurance claims made to the Medicare national social insurance program during the period 2006-2011 in addition to the quantity of ICD surgeries at the major cardiovascular referral centers known as cardiac care facilities (CCFs). The preeminent work in this domain is the Dartmouth Atlas in which the then Center for the Evaluative Clinical Sciences, Dartmouth Medical School [28] assigned hospitals to one of 306 health referral regions (HRRs) representing markets for tertiary medical care. The significant contributions made by the Dartmouth Atlas to health services research motivates our work but our methodology, which is based on network-graph topology, is a departure from the HRRs strict adherence to geographic proximity that continues to be pursued by some efforts to modernize the designations. [29] Related work in defining regions according to network topology for the purposes of measuring health care variation across regions are nonetheless fundamentally based on geography. [14, 17, 15] Note that geographic proximity could be an alternative or additional constraint to our network-based method but is not necessary. In particular, if the analysis health care services provided by telemedicine or remote monitoring is desired then geographic information may be ignored intentionally in the discovery of network communities.

We develop in the following a recursive backtracking procedure for greedy search of high modularity communities that are each constrained by the requirement to contain at least one of the so-called special vertexes, which in the present application are the network-graph representatives of the hospital network cardiac care facilities at which a minimum number of ICD surgeries were performed and apply the method to approximate the optimal health care community assignment for each hospital in the network.

The inclusion of constraints in community detection, particularly a constraint of the nature under present consideration, is pertinent across many domains and applications. The community structure in any social network consisting of entities possessing distinct classes or attributes, e.g., standard entities versus noteworthy entities, may be more accurately characterized by our methodology. For example, a researcher seeking to partition the individuals of social communication networks that are the result of correspondences among senior and junior members of an organization, e.g., the Enron corpus [19], may enforce through our methodology the constraint that each community of organization members in the network consist of at least one senior member. (A reason for doing this is that the researcher knows how the organization structures its workforce but does not know which senior employee is grouped with which junior employees). This example is retrospective in that latent groups are being rediscovered by the researcher. The same example could apply in a prospective form. For example, given a network of professional relationships among its employees of different ranks, a company could use the algorithm to form optimal groups with at least one senior employee per group. Analogous examples might also arise in peer mentoring situations – a teacher in a classroom might use a friendship network among the students to form workgroups such that each group consists of at least one student who is performing very well (on course for the top possible grade) who can help the other students understand the material. On the other hand, a researcher investigating the community structure in a trophic network [11] may desire a partition of species in which each community contains at least one member of the Canidae taxonomic family. We emphasize that applications of our methodology abound in networks consisting of heterogeneous entities.

We represent the nation-wide hospital network with the weighted, undirected graph $\mathcal{G}=(\mathbf{V},\mathbf{E})$ which designates one vertex $v\in \mathbf{V}$ for each hospital and a positively weighted edge $\{u,v,{w_{uv}}\}\in \mathbf{E}$ for each pair of vertexes $\{u,v\}\in {\mathbf{V}^{2}}$ that represents interacting entities in the network, e.g., patient-sharing hospitals, where ${\mathbf{V}^{2}}=\mathbf{V}\times \mathbf{V}$ is the set of vertex pairs. Note that if ${w_{uv}}=0$ then $\{u,v,0\}\notin \mathbf{E}$. The weight ${w_{uv}}$ of edge $\{v,u,{w_{uv}}\}$ reflects the quantity of shared patient visits recorded in Medicare claims data between the hospitals represented by vertexes $u,v\in \mathbf{V}$.

Suppose that the vertex subset ${\mathbf{V}^{\prime }}\subseteq \mathbf{V}$ contains the vertexes that are called special vertexes and represent network entities of a distinguishing nature. In the present scenario, ${\mathbf{V}^{\prime }}$ consists of the vertexes that represent the cardiac care facilities at which at least τ ICD surgeries, for some $\tau \ge 0$,2 were performed and define ${\mathbf{V}^{\prime }}={\mathbf{V}^{\prime }_{\tau }}$ accordingly. The information contained in the labeling of these special vertexes in ${\mathbf{V}^{\prime }_{\tau }}$ is information not encoded in the graph topology induced by the edge set E itself and, therefore, its integration into a standard community detection procedure must be made explicitly.

A useful mathematical representation of the weighted, undirected network-graph $\mathcal{G}$ in a variety of approaches to network science applications is the non-negative, symmetric $p\times p$ matrix W in which ${W_{ij}}={w_{{v_{i}}{v_{j}}}}$, for ${v_{i}},{v_{j}}\in {\mathbf{V}^{2}}$, i.e. $i,j\in \mathbb{Z}$ such that $1\le i\lt j\le p$ corresponding to some ordering of the vertexes ${v_{k}}\in \mathbf{V}$. The network modularity quality function where $\mathbf{s}\in {_{p}}={\{1,2,\dots ,p\}^{p}}$ is a community assignment vector consisting of at most p unique community membership labels, ${d_{j}}={\textstyle\sum _{i=1}^{p}}{W_{ij}}$ is the degree of vertex ${v_{j}}$, and $2m={\textstyle\sum _{j=1}^{p}}{d_{j}}$ is the total degree of the network-graph $\mathcal{G}$. Note that ${W_{ij}}$ is the observed (true) weight of the edge connecting vertexes ${v_{i}}$ and ${v_{j}}$ whereas ${E_{ij}}=\frac{{d_{i}}{d_{j}}}{2m}$ is the expected weight of the edge connecting these vertexes under randomization via the configuration model. [3]

The process of unconstrained modularity optimization involves the identifying of a partition $\{{\mathbf{V}_{1}},\dots ,{\mathbf{V}_{n}}\}\subseteq \mathbf{V}$ of size $n\le p$ corresponding to the community assignment vector s such that ${s_{k}}=r$ if ${v_{k}}\in {\mathbf{V}_{r}}$, for some $r\in \{1,2,\dots ,n\}$. For a partition of size $n\le p$, the Stirling number of the second kind [1] counts the number of unique community assignments. Note that if the ${p^{\prime }}=|{\mathbf{V}^{\prime }_{\tau }}|$ special vertexes are partitioned into n vertex sets then the number of feasible community assignments is with $n\le {p^{\prime }}$, which displays asymptotics similar to the number of unconstrained partitions. The special case when $n={p^{\prime }}$, for which $C(p,{p^{\prime }},n)={n^{p-n}}\in O({n^{p}})$, amounts to the unconstrained optimization of the modularity quality function over the set of non-special vertexes $v\in \mathbf{V}\setminus {\mathbf{V}^{\prime }}$ with each of the special vertexes $v\in {\mathbf{V}^{\prime }}$ each belonging to a unique community. This observation accordingly supports a bottom-up approach that is a part of our forthcoming procedure.

On the extreme end of the spectrum, if $\tau =\omega $, where ω equals the maximum number of ICD surgeries performed at any cardiac care facility in the network, then ${p^{\prime }}$ achieves its minimum tenable value. In particular, if the maximum number ω of ICD surgeries performed at any hospital in the network is unique then ${p^{\prime }}=1$. Moreover, if $\tau \gt \omega $ then no feasible solution to the constrained optimization problem exists. In general, ${p^{\prime }}$ is a non-increasing function of $\tau \ge 0$.

The community assignment vector that optimizes the modularity quality function where ${_{p}}={\{1,2,\dots ,p\}^{p}}$, is the community assignment vector which, on average, labels similarly vertexes that are more well-connected than expected to the same community among the $n\le p$ different communities.

Suppose that $R(\mathbf{s}|{\mathbf{V}^{\prime }_{\tau }})$ is the boolean function that returns true if the community assignment vector s satisfies the desired property that each community contain at least one special vertex ${v^{\prime }}\in {\mathbf{V}^{\prime }_{\tau }}$ and define the restricted community assignment vector that optimizes the modularity quality function where, in the present scenario, where $\mathcal{U}(\mathbf{s}|{\mathbf{V}^{\prime }_{\tau }})$ is, for example, the set of unique community labels for those vertexes $v\in {\mathbf{V}^{\prime }_{\tau }}$. Defined as such, $R(\mathbf{s}|{\mathbf{V}^{\prime }_{\tau }})$ returns $\text{True}$ when each community is constituted by at least one special vertex ${v^{\prime }}\in {\mathbf{V}^{\prime }_{\tau }}$.

We define a health care community assignment as the designation of hospitals to communities in a community assignment vector that optimizes the network modularity quality function $Q(\mathbf{s}|\mathbf{W})$ in Equation (2.1). The quantity $Q(\mathbf{s}|\mathbf{W})$ is proportional to the difference between the fraction of (weighted) edges within communities and the expected fraction if edges were randomly distributed according to the configuration model, i.e. edge randomization with degree distribution conserved, subject to the constraint that to each community belongs at least one special vertex, e.g. a vertex representing a cardiac care facility where at least $\tau \ge 0$ ICD implantations occurred. Specifically, we seek to maximize modularity while requiring that the number of cardiac care facilities per health care community is at least $\tau \ge 0$: a restriction we encode with Boolean variable $R(\mathbf{s}|{\mathbf{V}^{\prime }_{\tau }})$ indicating the feasibility of a community assignment.

Existing work in this realm considers incorporating additional information in the form of individual entity labels and pairwise constraints, i.e., that two vertexes must be labeled similarly or differently. [7] A restriction of the type we consider here has, to our knowledge, remained unstudied. In the context of the hospital network, the HRRs defined previously associated local health care markets to the tertiary care facilities where the plurality of the residents were referred for major cardiac procedures. An HRR is a reflection of its regional health care market and, because necessarily within each is a hospital specialized in cardiac surgery, a comparison across regions is facilitated. In an effort to standardize cardiac care among health care subnetworks, we present an initial undertaking towards developing a paradigm of constrained community detection. In particular, because constraint satisfaction in optimization problems provides context, the health care communities identified by our constrained community detection approach have real-world utility.

The organization of this article is as follows. In Section 3, we mathematically formulate the constrained optimization problem. In Section 4, we define a greedy, recursively-backtracking procedure for identifying high-quality, constrained communities. Utilizing our procedure in Section 5, we estimate the health care communities of the nation-wide hospital network and illustrate in Section 6 our method on a network-graph of well-known community structure.

Suppose that $\mathcal{G}=(\mathbf{V},\mathbf{E})$, with $p=|\mathbf{V}|$, is a network-graph and that ${\mathbf{V}^{\prime }}\subset \mathbf{V}$ is a subset of vertexes. Moreover, suppose that $\mathbf{s}\in {\{1,2,\dots ,n\}^{p}}$, for some $n\in {\mathbb{Z}_{+}}$, is a vector of labels such at each label is applied at least once. We define the collection of $\left(\genfrac{}{}{0.0pt}{}{p}{2}\right)$ binary variables ${x_{ij}}=1\{{s_{i}}={s_{j}}\}$, for $i,j\in \{1,2,\dots ,p\}$, along with the corresponding collection of edge values ${b_{ij}}={w_{ij}}-{d_{i}}{d_{j}}/2m$ and the vertex values ${u_{k}}=1\{{v_{k}}\in {\mathbf{V}^{\prime }}\}$, for $k\in \{1,2,\dots ,p\}$. Let $\mathbf{B}=[{b_{ij}}]$.

The constrained optimization problem under consideration is described as follows. This discrete optimization problem is, in general, at least as difficult as the corresponding, unconstrained NP-hard decision problem which seeks an answer to whether or not a partition of the vertex set exists with a modularity quality function value of at least some minimum value. In particular, the boundary value $\eta =0$ corresponds to the special case of an unconstrained optimization problem in the description above.

A greedy unconstrained method may readily agglomerate communities toward achieving a local maximum of the network modularity function $f(\mathbf{s},\mathbf{B})$. On the other hand, satisfying the constraint requires the consideration of many permutations of such mergers and any recursive backtracking procedure for identifying high-quality, with respect to the network modularity function, has potentially intractably-many paths in the search space to traverse. Note that if the graph $\mathcal{G}$ happens to be unweighted then the adjacency matrix W is binary and, up to quantities involving a $2m$ denominator, Equations (4.1) and (4.2) involve integers on similar scales. Accordingly, the various sequences of community assignment label updates through the recursive backtracking procedure are reduced. In the application of our method toward the defining of health care communities in the nationwide social network of cardiac-related Medicare referrals between hospitals we make this binary transformation $\mathbf{A}=1\{g(\mathbf{W})\gt \zeta \}$, where $g:{\mathbb{R}^{p\times p}}\mapsto {\mathbb{R}^{p\times p}}$ is defined in Section 5.1, the indicator function $1\{P\}=1$ if proposition P is true and $1\{P\}=0$ otherwise, and $\zeta \gt 0$ is determined from the data, to facilitate tractable computations.

We begin this section by presenting the existing work on unconstrained network modularity quality function optimization upon which our generalized procedure, which seeks to maximize the network modularity quality function over the space of feasible community assignment vectors, is based. Subsequently we outline our method for constrained optimization of the network modularity quality function and provide pseudocode for pertinent procedures.

The nationwide hospital network we consider consists of $p=4734$ hospitals, as depicted in Figure 1a, among which 1388 (29.3%) are CCFs hosting at least $\tau =1$ ICD surgeries and correspond to special vertexes in ${v^{\prime }}\in {\mathbf{V}^{\prime }}$. We begin this analysis by processing the data.

In the following, we consider each k-tuple, for $k\in \{2,3,4,5,6\}$ as the set of special vertexes ${\mathbf{V}^{\prime }}$, among the set of $p=34$ vertexes in the unweighted Zachary Karate Club social network-graph. [32] We apply Procedure 1 to determine, for each of these tuples of special vertexes, to record the modularity of each constraint-satisfying community assignment vector returned by Procedure 1, the initial position selected by the procedure, and the relative length of the computational time for the procedure to halt. The results of the simulation study are presented in Figure 4.

The quantity $C(p,{p^{\prime }},n)$ in Equation (2.2) counts the number of feasible community assignments for a given $p=|\mathbf{V}|$, ${p^{\prime }}=|{\mathbf{V}^{\prime }}|$, and number of communities n. While this number reflects the number of community assignments necessary to brute-force check and, therefore, guarantee that the optimum defined in Equation (2.4) has been obtained, our greedy procedure is guided by the network topology and halts in many fewer iterations. Since the number of communities n is automatically chosen by the procedure, the computational time required for our procedure to halt is a function of (i) the number ${p^{\prime }}$ of special vertexes and (ii) the distribution of the special vertexes within the network-graph. For example, if ${n_{opt}}$ is the number of unique labels, i.e. communities, represented in ${\mathbf{s}_{opt}}$ in Equation (2.3) and ${p^{\prime }}\lt {n_{opt}}$ or if ${p^{\prime }}\ge {n_{opt}}$ but the special vertexes are frequently labeled similarly in ${\mathbf{s}_{opt}}$ then some work is necessary to compute ${\mathbf{s}_{R}}$ of in Equation (2.4).

Our method for identifying network communities optimizes a quality function while adhering to constraints. The results in this paper establish the utility of penalized optimization in community detection. Our method is versatile and amenable to many types of constraints on the composition of communities. We note that our procedure is valid for any constraint which is an increasing function of the variable of interest, e.g., number of CCF hospitals belonging to a community, number of cardiac surgeons, quantity of cases involving improper medical procedures, etc. The key requirement of our constrained optimization procedure is that the merging of two communities must not be the basis for the resulting community to be in violation of the constraint. A constraint that imposes a maximum ICD volume is, for example, not of this type.

There exists a disconnect between network science and health services research due in part to the incongruence between mathematical elegance and real-world constraints. We have provided an illustration of the application of both a pure (unconstrained) method and one with constraints. We solved the practical problem of partitioning a network of hospitals with the constraint that the number of ICD surgical procedures that have taken place at at least one hospital belonging to each community exceeds some threshold. Though our method advances both the community detection and heath services literature, it is not complete from the perspective of a health care policy maker since many real-world constraints remain to be incorporated. Another type of constraint is, for example, given the geographic locations of hospitals, a requirement that communities not exceed a defined geographic maximum diameter or that they satisfy a geographic congruity constraint. On the other hand, one shouldn’t necessarily seek to impose geographic contiguity constraints, for example, if one is interested in analyzing the effect of telemedicine or remote monitoring, for which the organization of health care does not need to conform as much to geography. This article is an initial exploration of a line of thinking that we anticipate will substantially advance the practical utility of community detection.

In terms of health policy, our future research involving an outcomes-based analysis of the communities discovered by our method, as constrained here by minimum ICD surgery volume and subsequently by other factors, will lead to enhanced acuity and potentially greater statistical power for studying variations in health care markets (based on patient referral patterns). Through standardizing the composition of the HCCs, our method provides the tools for such comparisons to be made meaningfully.