Spatial regression and spillover effects in cluster randomized trials with count outcomes

Funding information MedicalResearchCouncil,Grant/Award Numbers:G7508177,MR/K012126/1; EuropeanandDevelopingCountriesClinical Trials Partnership,Grant/AwardNumber: EDCTP2;Department for International Development,Grant/AwardNumber: MR/K012126/1 Abstract This paper describes methodology for analyzing data from cluster randomized trials with count outcomes, taking indirect effects as well spatial effects into account. Indirect effects are modeled using a novel application of a measure of depth within the intervention arm. Both direct and indirect effects can be estimated accurately even when the proposed model is misspecified. We use spatial regressionmodels with Gaussian random effects, where the individual outcomes have distributions overdispersedwith respect to the Poisson, and the corresponding direct and indirect effects have a marginal interpretation. To avoid spatial confounding, we use orthogonal regression, in which random effects represent spatial dependence using a homoscedastic and dimensionally reducedmodification of the intrinsic conditional autoregression model. We illustrate the methodology using spatial data from a pair-matched cluster randomized trial against the dengue mosquito vector Aedes aegypti, done in Trujillo, Venezuela.


INTRODUCTION
A cluster randomized trial (CRT) is an experiment in which groups, rather than individuals, are randomized (Hayes and Moulton, 2009). The clusters are usually defined geographically (eg, houses or communities) or institutionally (eg, health facilities or schools). One reason for randomizing by cluster is to achieve a mass or indirect effect (Halloran and Struchiner, 2009), which we take to mean an effect exerted by the intervention on individuals who do not receive it. Another reason is to reduce contamination, that is, some individuals in the control arm receiving an intervention intended for the other arm (Hayes and Moulton, 2009): this is not addressed in the current paper.
We identify indirect effects resulting from proximity between clusters, which we call spillover indirect effects. We in the current sense, is important for some classes of intervention, and is a specific objective of some trials (Hawley et al., 2003;McCall, 2011).
We also discuss another type of spillover effect. Independence of the clusters may be violated if they are not sufficiently separated, analogous to agricultural interplot interference (Besag and Kempton, 1986). Ideally, betweencluster independence would be achieved by having sufficiently distant clusters. For example, McCann et al. (2018) propose a selection method for candidate clusters that minimizes cross-cluster effects (they use the term contamination for such effects, a definition that is more general than the one stated above). However, it is not always feasible to guarantee that the clusters are independent for practical purposes, yet frequently no allowance is made for this in analysis (Jarvis et al., 2017).
We discuss two types of possible dependence: intracluster correlation (Hayes and Moulton, 2009) and dependence between individual outcomes in different clusters. We call this last type spillover dependence. This dependence is assumed to exist prior to the intervention and may arise from: (a) the spatial configuration of the individuals in the trial, (b) contact patterns of infectious diseases and (c) individuals belonging to more than one cluster either sequentially or simultaneously (Silcocks and Kendrick, 2010).
For two-arm CRTs with a count outcome, we propose spatial regression methods to identify, and adjust for, both types of spillover effect, while incorporating the clusterrandomized design. Spillover indirect effects are addressed via an individual-level covariate that measures the degree of surroundedness of individuals in either arm to others in the intervention arm. The estimates of the intervention effects, and of the spillover indirect effects, have simple interpretations in terms of marginal expectations and they are well defined even when the proposed model does not hold. Spillover dependence is addressed via spatially correlated random effects.

The proposed individual-level model
When individual level count data are available, generalized linear Poisson models with random effects are commonly used to analyze CRTs, with one random effect per cluster (Hayes and Moulton, 2009). Suppose there are clusters and a total of individuals in the trial. The count outcomes are denoted by 1 , 2 , … , and the corresponding denominators by 1 , … , , while = ( 1 , … , ) is the vector of cluster random effects ( denotes the transpose of ). In this paper, we assume that each outcome is associated with a single geographical point location and that such information is available. We use = ( 1 , 2 , … , ), a vector of spatially correlated random effects, to model geographical dependence. Conditionally on the random effects, 1 , … , are assumed independent with Poisson distribution and conditional means given by ( | , ) = exp{ + + The binary variable indicates whether individual is in a control or intervention cluster ( = 0, 1, respectively). The variable and its interaction with are introduced to model spillover indirect effects (Section 2.3). The standard regression model for CRTs with count outcomes (Clark and Bachmann, 2010) corresponds to the case where there are no spillover indirect effects ( = = 0) and also no geographical dependence ( 2 = 0 so that ≡ 0), namely where the 0 subscripts are introduced to distinguish the parameters from those in (1). We will call Equation (1) the extended model. We use ( , ) to denote the multivariate normal distribution of dimension with zero mean vector and covariance matrix , and denotes the identity matrix of dimension × . We assume the cluster random effects vector has the form , where ∼ ( , 2 ) and the cluster design matrix contains cluster-level covariates, where the th covariate is a vector with 1s in entries corresponding to individuals in cluster , and 0s elsewhere. To model the geographical dependence, we restrict the spatial random effects to have linear form, namely = , where is a known matrix with columns that takes into account the spatial structure and ∼ ( , 2 ). The linear form assumption, and the choice of (the number of columns of ), are discussed in Sections 2.4.3 and 2.4.4.
Let denote the fixed effects design matrix of the extended model (1). Then the linear predictor can be written in vector form as + + , where = ( , , , ) . We restrict the spatial design matrix to satisfy the condition = so its columns are orthogonal to those of . The orthogonality between fixed and random effects design matrices is required to avoid so-called spatial confounding that is discussed in Section 2.4.1.

Marginal interpretation
By integrating out all the random effects, the extended model (1) can also be thought of as a marginal probability model, known as the multivariate Poisson-lognormal distribution (Aitchison and Ho, 1989). There is no closed-form expression for the density but the marginal means can be expressed as (3) where we define ∶= Var( + ) = 2 + 2 and is the th row of . The standard model can also be seen as a marginal model but where = Var( ) = 2 does not depend on . The marginal variance of the extended model is Var( ) = ( )[1 + ( ){exp ( ) − 1}] which shows that the marginal distribution of is overdispersed with respect to the Poisson.

Intervention effects
In this section, we identify different ways to quantify the effect of the intervention in the population. Their use will depend on the aims of the particular cluster randomized trial. We emphasize that all the intervention effects defined in the following sections (eg, 2.2.1, 2.2.2, and 2.3) are meaningful irrespective of any underlying assumed model. This is because they are defined in terms of marginal expectations of the outcome variable . In particular, these intervention effects are meaningful even when the extended model does not hold.

Pairwise intervention effect
We define the individual-level pairwise intervention effect as where and correspond to individuals in the intervention and control arm, respectively. In particular, for the standard model (2) we have that Pint = 0 for all , . However, for the extended model, Pint = + − + ( − )∕2, which varies over ( , ). Specifically, there is dependence on ( , ) as well as on ( , ).

Total intervention effect
Let 1 and 0 be the set of indices of individuals in the intervention and control arms, respectively, and let + = ∑ ∈ for = 0, 1. The total intervention effect is defined as that is, the logarithm of the ratio between the expected raw rates in the intervention and control arm. There are alternative choices for a total intervention effect (see Chapter 10 of Hayes and Moulton, 2009 and the discussion section in this paper) and we do not evaluate them in this paper. For the extended model (1) the total intervention effect Tint can be written as which from now on we denote by , and where ( , , 2 ) = log which does not depend on 2 as it cancels in (7). Based on its factorization in Equation (6), we can interpret as the sum of two terms: , the direct contribution of each intervention and individual; and , the contribution of all the indirect effects. In Section 2.3.4, we show that and keep the same interpretations outside the extended model. Under the standard model ( , , 2 ) ≡ 0 and therefore Tint = 0 is the total intervention effect. This highlights the fact that, perhaps contrary to intuition, the parameter , when the extended model holds, has a different interpretation to that of 0 when the standard model holds, both being the multipliers of the intervention variable in each case. As opposed to the pairwise intervention effect Pint that in general depends of the chosen pair ( , ), we emphasize the statistical inference regarding the total intervention effect Tint, which is a single quantity and is meaningful even under misspecification of either the standard or extended model. Estimates of Tint under the standard model are simply estimates of 0 while, under the extended model, they are estimates of , constructed using expressions (6) and (7). We describe our estimation procedures in Section 3.

The constant variance property
A special case of the total and pairwise intervention effects in the extended model arises when the variances Var( ) = 2 do not depend on . In this situation, both the total intervention effect parameter and all the pairwise intervention effects depend only on fixed effects , , and and not on any variance component 2 or 2 . This motivates the following definition: Definition: The constant variance (CV) property A model of the form (1) has the constant variance property when the spatial variance Var( ) does not depend on .
A sufficient condition for this property is that the rows of the spatial design matrix have unit norm, that is, = 1 for all . Henceforth, we further restrict our proposed extended model (1) so that the rows of spatial design matrix have unit norm and therefore the parameter in (7) only depends on and and simplifies to In Section 2.4.2, we describe an algorithm to make the unit norm restriction on the rows of consistent with the orthogonality condition on its columns.

Spillover indirect effects
CRTs usually include a small number of clusters (perhaps around 20 or 30) and the allocation may be geographically imbalanced in the sense of special spatial configurations. Moreover, some trials have clusters that abut each other, for example, divided by only a narrow street. We are interested in quantifying the indirect effect that such configurations may exert on trial outcomes. Here the question arises of whether intervention locations exert any indirect effect only on control locations, or also on fellow intervention locations. For generality, we allow the effects of intervention locations to differ according to the arm of the nearby locations, and hence define:

Contralateral spillover indirect effect: An intervention effect exerted on individuals in the control arm via nearby individuals in the intervention arm;
Ipsilateral spillover indirect effect: An intervention effect exerted on individuals in the intervention arm, accentuated via nearby individuals that are also in the intervention arm.
To represent these types of indirect effects in the extended model (1), we need a measure of the degree of closeness or surroundedness that each individual spatial location has, with respect to intervention individuals. Two such measures are described next.

Measures of surroundedness
As a measure of surroundedness, we could use the linear distance from an individual (in either arm) to the closest of those in the intervention arm. However, this does not capture a possible mass effect exerted by, for example, multiple individuals in the intervention arm surrounding one in the control arm. Hence we propose that any surroundedness variable should be defined as the number of intervention individuals who, in some sense, are around the individual . For example, we can use the number of individuals within a specified radius who are in the intervention arm (henceforth the disc measure). We also propose depth (Tukey, 1975) as a measure that does not require a radius to be specified. In general, depth measures the spatial centrality of one point relative to a given set of others. To our knowledge, statistical depth has only been used as a measure of multivariate location, so its use in the context of spatial indirect effects is novel. Here we will assume only two dimensions are relevant and work in the plane, using Tukey's half space depth: rotating a doubly infinite line through a spatial location of interest, and monitoring the number of spatial locations on each side of it, Tukey's half space depth of the index location is the smallest count occurring on either side as the line is rotated through 180 • . This can be computed using the algorithm of Rousseauw and Ruts (1996). For an individual in the control arm, depth is computed with respect to the set of all intervention individuals. For an individual in the intervention arm, the depth is computed with respect to the set of intervention individuals excluding itself. Either of these proposed measures of surroundedness is included in the model as a fixed effect. We elaborate further on these measures in the context of specific geographical configurations used in the simulations in Section 4.

Pairwise spillover indirect effects
The pairwise spillover indirect effect is exerted by surrounding individuals, acting on the expected individual marginal rates of the outcomes. Although such an effect may be exerted across arms (contralateral), the gradient is defined with reference to individuals in the same arm. So, letting and be such individuals ( = ) for which − = 1, then the pairwise spillover indirect effect is defined as For clarity, we define Sind 0 as the pairwise spillover indirect effect for pairs of individuals in the control arm and Sind 1 as the pairwise spillover indirect effect for pairs of individuals in the intervention arm. Under the extended model, and assuming the CV property, we can write Sind 0 = and Sind 1 = so they do not depend on the particular pair of individuals selected, but only on whether the individuals are both in the control or both in the intervention arm. For each extra surrounding individual in the intervention arm, individual 's expected rate changes by a factor of exp( ) if is in the control arm, and by exp( ) if in the intervention arm. In this way, then represents the pairwise contralateral spillover indirect effect, and the pairwise ipsilateral spillover indirect effect.

Total spillover indirect effects
Let 1 and 0 be the number of individuals in the intervention and control arms, respectively. To define a total contralateral indirect effect, we first note that we can assume that the values of the surroundedness variable for individuals in the control arm satisfy 1 ≤ 2 ≤ ⋯ ≤ 0 . This is always possible via relabeling of the indices in 0 . The total contralateral indirect effect (due to the surroundedness ) can then be defined as that is, the average of all pairwise contralateral indirect effects. The ordering restriction on the surroundedness variable is imposed to be able to interpret this effect as the average effect that increasing the surroundedness variable has on individual rates. The total ipsilateral indirect effect Tind 1 can be defined in a similar fashion. In the extended model, the total contralateral and ipisilateral indirect equal Δ 0 andΔ 1 , respectively, wherē so they are simple multiples of and , respectively, with corresponding factors equal to the average pairwise difference of the surroundedness variable in the control and intervention arm, respectively. In this way, we can interpret total contralateral and ipsilateral effects geometrically where, for example, the effect of 10 extra surrounding intervention individuals reduce the individual expected response by a factor of exp(10 ) for individuals in the control arm.

Total residual indirect effect
In this section, we give interpretations to the parameters , , and in the extended model (1) that are valid even when this model does not hold. Any proposed measure of surroundedness will be a discrete variable that allows individuals to be defined as isolated (from the intervention arm) if = 0. Let 1 and 0 denote the set of indices corresponding to isolated individuals in the intervention and control arm, respectively. We now define the total intervention effect for isolated individuals: where + = ∑ ∈ for = 0, 1, that is, the logarithm of the ratio between the expected rate of isolated individuals in the intervention arm and the expected rate of isolated individuals in the control arm. In the extended model, the total intervention effect for isolated individuals equals . We can now give an interpretation to the parameter defined in (8) that is well defined outside the extended model. We define the total residual indirect effect Tred simply as the difference between the total intervention effect Tint and the total intervention effect for isolated individuals Tiso, that is Evidently, in the extended model, the total residual indirect effect Tred is simply given by , since = − . Expression (13) reveals the intuitive fact that the difference between a total effect (Tint) and a pure effect (Tiso) is a measure of indirect effect. In practice, however, it is possible that there are no isolated individuals in either the intervention or the control arm, that is, either 1 = ∅ or 0 = ∅. In that case, as there is no information in the data to estimate these effects without an assumed model, we say that the Tiso and the Tred are not estimable. Finally, we define the total average expected count for control isolated individuals as This quantity is equal to exp( ) in the extended model, and exp( 0 ) in the standard model.

Spatial modeling
The main reason to incorporate spatial dependence in our extended model is that the randomization units (clusters) of many trials are defined geographically. In the CRT setting, the main interest is in the estimation of direct and indirect effects of the intervention, with spatial dependence being of interest insofar as it affects such estimation. Technically, the spatial random effects vector is used to implement a smoother, with the shrinkage estimates of the fixed effects ( , , , ) being interpreted accordingly (see Hodges, 2013, Chapter 10). For this reason, we focus the specification of the spatial dependence within the constrained setting that we have adopted. In this setting, = and is a matrix to be specified that should satisfy the CV property and should be orthogonal to the fixed effect design matrix in order to avoid the so-called spatial confounding that is described next.

Spatial confounding
Adding spatially correlated random effects such as in the extended model (1) can induce bias in the estimations, and artificially inflate the variances, of the fixed effects of interest ( , , , and ) (Reich et al., 2006;Hodges and Reich, 2010) directly or through the estimation of any of the total effects defined in Sections 2.2.2 and 2.3.3. One reasonable procedure is to retain only random effects that are orthogonal to the fixed effects matrix . Specifically, let = ( 1 ∕ 1 , … , ∕ ) be the random vector of individual rates, then the extended model (1) can be written in vector form as log{ ( | , )} = + + , where the × 4 design matrix is assumed to have rank 4, and is a proposed spatial design × matrix, which we assume has rank < − 4. Spatial confounding happens when the spatial random effects variance 2 is large with respect to the residual variation of the model (Hodges and Reich, 2010) and, at the same time, the eigenvectors of (the covariance matrix of ) associated with the largest eigenvalues, are strongly correlated with .
Following Hughes and Haran, 2013, there exist orthogonal matrices and of dimension × 4 and × ( − 4), respectively, such that = + . Then we can write + + = + 1 + + 2 , where 1 = and 2 = . The spatial confounding comes from the fact that and span the same column space. We follow the procedure of Reich et al. (2006) who suggests deleting the random effect term 1 altogether from the model in order to avoid spatial confounding. The surviving term is 2 = ∼ ( , ), where the covariance matrix = is rank deficient, with rank equal to min( , − 4) = . We let be the × matrix, with the eigenvectors of associated with nonzero eigenvalues, and normalized so that each eigenvector has norm equal to its corresponding eigenvalue. This implies that the constructed has the same dimensions as the initial , and 2 has the same distribution as , where ∼ ( , 2 ). It is easy to show that we can take the eigenvectors of to be orthogonal to , implying that the constructed satisfies the orthogonality condition = .
We note that the presence of an intercept term in the extended model means that the spatial random effect vector = , constructed above, satisfies the sum-to-zero condition: ∑ =1 = = = 0. The orthogonalization implicitly assumes that most of the variability in the outcomes is due to the covariates and not the random effects, so that the latter merely introduce spatial correlation between the outcomes. We should therefore expect the point estimates of the fixed effects under the orthogonalized model to be similar to those under the simpler model with no random effects at all (Reich et al., 2006).

Construction of orthogonalized models with the CV property
There is a compromise between the column orthogonality condition on and the unit-norm condition on its rows (the CV property, Section 2.2.3). We provide a simple algorithm to transform an arbitrary × matrix into a matrix of the same dimensions that satisfies both conditions and which is close to the original . First, a spatial random effects vector of the form = 0 has the CV property if the rows of 0 are of unit length so that the diagonal of the covariance matrix Var( ) = 2 0 0 has all entries equal to 1. Transforming 0 to 1 = {diag( 0 0 )} −1∕2 0 ensures that the resulting vector = 1 has the CV property. Now, for the orthogonalization, we proceed as in Section 2.4.1 by constructing 2 with the normalized eigenvectors of the matrix 1 1 . The algorithm consists in repeating, in turn, the transformations to the CV property and to orthogonality, and stops when the maximum absolute angle difference (to 90 • ) is within a predefined small tolerance. This is an iterative projection algorithm (Bauschke and Borwein, 1996) that targets the intersection between two sets so that if the two sets have nonzero intersection, the iteration converges linearly with the total number of individuals .

Choice of the spatial design matrix
There are random effect spatial models available that already have the required linear structure = , such as the one proposed by Silcocks and Kendrick (2010) that addresses the situation of membership of multiple clusters. However, here we are largely motivated by trials that randomize geographically (eg, houses based on their location) in which single cluster membership can be assumed. Other powerful and widely available spatial modeling strategies such as geostatistical modeling (Diggle and Ribeiro, 2007) for georeferenced data and intrinsic conditional autoregressive (ICAR) models for aereal data, need to be adjusted to satisfy our random effects specific structure.
ICAR models (Reich et al., 2006) use only the spatial neighborhood information in an × adjacency matrix . If the corresponding set of neighbors forms one connected component, then we can write = 1 √ + icar , where ∼ −1 ( , 2 −1 ) and = ∑ =1 has a flat improper distribution. The matrix icar is of dimension × ( − 1) and its columns are the eigenvectors of = diag( ) − that are associated with positive eigenvalues and also normalized so that their length equals the reciprocal of its corresponding eigenvalue. All the columns in icar are orthogonal to , so imposing the sum-to-zero constraint ∑ =1 = 0 yields the proper model with linear form = icar . This is equivalent to removing the spatial confounding with respect to only the intercept fixed effect .
Under the sum-to-zero constraint, it is easy to show that the variances of = icar are not constant as a function of so that the ICAR specification above does not have the CV property. After enforcing this property using the algorithm in Section 2.4.2, the resulting random effect is no longer ICAR. Geostatistical models, conversely, naturally satisfy the CV property and can be written in the form = geo (Lindgren et al., 2011), but the orthogonality condition needs to be enforced. Furthermore, the covariance matrix in geostatistical models depends on 2 as well as a scale and a smoothness parameters so, unless those extra parameters are known, the algorithm in Section 2.4.2 cannot be used and further modification is needed (see Hodges and Reich, 2010, for a related discussion).

Dimension reduction in the spatial model
We now discuss the choice of the number of columns of the spatial design matrix. The number of columns in a selected can be substantially reduced by discarding columns in the span of that give rise to negative dependence when there is no interest in modeling this (Hughes and Haran, 2013), as is the case in many CRT applications. The choice is guided by the spectrum of the modified Moran's operator ⟂ ⟂ (Hughes and Haran, 2013), where is the adjacency matrix of the neighborhood structure and ⟂ is the projector matrix onto the space orthogonal to the column span of . Moran's operator ⟂ ⟂ appears in the numerator of a generalized form of Moran's statistic : ( ) ∝ ( ⟂ ⟂ )∕( ⟂ ) . The standardized spectrum of Moran's operator comprises all possible values of ( ), and the eigenvectors comprise all possible mutually distinct patterns of spatial dependence. Positive (negative) eigenvalues correspond to varying degrees of positive (negative) spatial dependence. If positive spatial correlation is of interest, then we can further restrict the spatial design matrix to have the additional property that the diagonal elements of ( ⟂ ⟂ ) are all positive. This can lead to a substantial reduction on the original number of columns (Hughes and Haran, 2013). In this way, an initial matrix can be transformed to a matrix with the desired three properties, using a modified version of the algorithm in Section 2.4.2 where the orthogonalization steps are followed by column selection where only columns that give rise to positive spatial dependence are selected. This last step preserves orthogonality.
In summary, we propose that, in the absence of a specific spatial design, one should start with ICAR and transform it, using the iterative procedure described above, into a column-reduced matrix that satisfies both the CV and orthogonality properties. As shown by Hughes and Haran (2013), the resulting models are attractive for spatial modeling, for example, in retaining only positive spatial associations and also with smooth patterns of spatial variation at various scales. As we assume point spatial data are available for each individual, we define the neigborhood structure from the corresponding Voronoi tessellation (Bivand et al., 2008). This partitions the continuous plane into a set of regions called tiles where each tile is a convex polygon, defined as the subset of the plane for which the given location is closer than any other. We propose to use queen type neighbors instead of rook type (Bivand et al., 2008) as the former are more consistent with our definition of surroundedness.

STATISTICAL INFERENCE
To fit the extended model (1) to a given dataset, we use Bayesian inference for generalized linear mixed models as described in Fong et al. (2010) and, more specifically, integrated nested Laplace approximations (INLA; Rue et al., 2009), which is an efficient alternative to MCMC in our our situation of only two hyperparameters ( 2 and 2 ). For the fixed effects, we assign a flat improper prior for the intercept and independent 1 (0, 0.001) priors for , , and . For the random effect standard deviations and , we use independent exponential priors. These are weakly informative priors (called penalized complexity priors, Simpson et al., 2017) that favor the standard model in order to avoid fitting artificial spatial and cluster dependence. This is in contrast to commonly used Gamma priors for the precisions (1∕ 2 and 1∕ 2 ) that do not allow very small values of the standard deviations and therefore can cause overfitting. Following Simpson et al. (2017), the penalty is given by the magnitude of the marginal standard deviation of the random effect, that is, after integrating out the prior.
INLA naturally produces marginal posterior densities for the fixed effects, which, in turn, implies marginal posteriors for the total contralateral and ipsilateral indirect effects that depend only on and , respectively. To compute marginal posteriors of the total intervention effect and the total residual indirect effect , which are nonlinear functions of , , and , we use the Monte Carlo sampling procedure implemented in the INLA package in R (The R-INLA Project, 2019). This generates independent samples from the joint posterior of the fixed effects (based on a Gaussian approximation) so that we can apply Equations (6) and (8) to obtain samples from the marginal posteriors of and . The approximation is particularly good in our case as neither nor depends functionally on the hyperparameters 2 or 2 .

SIMULATION STUDY
We performed a simulation study in order to evaluate the operating characteristics of our method. We performed six sets of simulations: one for each measure of surroundedness described in Section 2.3.1 (disc and depth, each included as a fixed effect) and one for each of three spatial configurations. In line with similar simulation studies (Paciorek, 2010;Hughes and Haran, 2013), we assume the spatial configuration of the trial is a unit square grid in which each square represents an individual and groups of squares represent clusters. Moreover, we use square clusters so that the neighborhood structure is obvious. In order to prevent artifactual behavior of the depth being induced by the regular grid layout of the square centroids (Rousseauw and Ruts, 1996), we perturbed the individual centroids slightly. We use a completely balanced design with 18 clusters in each arm and with nine individuals per cluster giving a total of 324 individuals. The three cluster configurations shown in Figure 1 were chosen to illustrate different patterns of surroundedness, rather than being "typical" spatial randomizations (although, as discussed in Section 6, this notion is not trivial to define).
With the aim of assessing the robustness of the proposed method, the generating model was not the extended model. However, an intervention effect and indirect effects were present (details are in the Supporting Information section). Within each set, two models were fit: the full extended model (1), and the standard model (2). The latter was introduced to measure the influence of ignoring spillover indirect effects when they are present in the data. In this situation, as discussed in Section 2.2.2, estimates of the parameter 0 in the standard model are actually estimating the total intervention effect Tint. Results for the disc measure of surroundedness are shown in Figure 2. For the three configurations, in terms of the posterior median, the method was able to correctly estimate the target parameters (thick vertical line): the total intervention effect Tint (when using the standard or extended model) and the total contralateral and ipisilateral indirect effects (Tind 0 and Tind 1 , respectively). However, for the Ring configuration, there seem to be a small bias in estimating Tint for the extended model that is corrected by using depth instead of the disc measure of surroundedness (corresponding histograms are in the Supporting Information section). This reflects the clear differences in values for both measures of surroundedness as shown in Figure 1. The only small bias that remains, when using both measures of surroundedness, is when estimating the total contralateral indirect effect Tind 0 . Other differences not shown by the histograms are structural and are related to estimating other parameters such as the total residual indirect effect (Tred) and the total intervention effect for isolated individuals (Tiso). When using the disc measure, there were no isolated individuals in the intervention arm in the three configurations considered. As described in Section 2.3.4, this means that Tred and Tiso are not estimable and therefore the corresponding estimates of the parameters and are still well defined but are estimating a different population parameter. This was not the case for the depth, as some individuals are defined as isolated under this measure.
In further simulations, we allowed the cluster sizes to vary, in terms of both numbers of individuals and spatial area. The conclusions remained the same, with no indication that this variation induced bias in the results.

APPLICATION: TRUJILLO TRIAL
In this section, we apply the above methods to a CRT of interventions against Aedes mosquitoes, which are vectors of dengue and Zika viruses. In Trujillo, Venezuela, Kroeger et al. (2006) randomized pairs of clusters of households to a package of window curtains and water container covers treated with insecticide, or to control (no intervention). A total of 18 clusters were defined geographically and pairmatched in a semi-quantitative manner, based primarily on baseline indices of infestation of water containers by F I G U R E 1 Three configurations used in the simulation study (leftmost column of three panels), and the depth and disc measures of surroundedness (middle and right columns, respectively). In the left column, the first configuration (top), which we call Ring, the control clusters (white) are almost completely surrounded by the intervention ones (gray). In the second configuration (middle), called Crater, control clusters are again contiguous, but tending to lie on one side of the space, with the interventions clusters on the other. In the third configuration (bottom), called Chessboard, control and intervention clusters alternate regularly over the study region. The middle column shows values of the depth measure of surroundedness. Darker values show higher values. The right column shows values of the disc measure, that is, number of individuals in the intervention arm who lie within a radius of each individual. Here the radius is set at 0.12, the study area having unit side. Again, darker values show higher values mosquito larvae and pupae, but also on residential conditions such as the size and density of houses. The pairing, together with the randomization within each pair, agglomerated the clusters into two contiguous regions per arm, with one control region being interdigitated with one intervention region. Point geo-reference information was available for each house site (latitude and longitude) and this was used to construct the corresponding Voronoi tesselation ( Figure 3). The primary outcome was the number of containers per house that were positive for immature stages of the mosquito. The house in this trial corresponds to the more general term "individual" as used in the current paper. When this outcome is aggregated, for example, to the clus-ter level, and expressed per 100 houses, this is called the Breteau index. The flight range of Aedes mosquitoes is in the hundreds of meters (Honório et al., 2003), which is the same order of magnitude as the size of the area of the current trial. Daily human movements within urban areas are likely to be larger, and relevant for dengue transmission (Stoddard et al., 2013). Hence people may be infected outside the study area, but such outcomes were not endpoints for the current trial. The final entomological survey was done 12 months after the interventions were implemented. The original report analyzed the cluster-level changes in Breteau index from baseline to the final survey, and found no significant difference between control and intervention arms (Kroeger et al., 2006).

F I G U R E 2
Simulation results. Histograms of posterior medians from 1000 independent replicates. Top, middle, and bottom rows correspond to Ring, Crater, and Chessboard configurations (see Figure 1), respectively, and all for disc measure of surroundedness ( = 0.12).

Thick vertical lines indicate population (target) parameters values
In a post hoc analysis, the original report found some evidence of a contralateral spillover indirect effect, in that initially positive houses (ie, with immature mosquito stages) in the control arm within 50 m of the nearest intervention house were 3.5 times as likely to be free of vectors than more distant initially positive control houses 1 month after the intervention, although this difference was not statistically significant (Kroeger et al., 2006).
The current re-analysis again uses the final survey data that included data from 702 house sites. Another independent random effect was included in the extended model (1) to account for the cluster pairing. For prior sensitivity, we tried three different priors for the standard deviation of the spatial random effects . The priors were set such that the corresponding marginal standard deviations were 0.3, 0.05, and 0.01, which corresponds to weak, medium, and strong penalties to the spatial model as compared to the model with no spatial random effects. For the cluster and pairing random effects, we only used (independent) priors with weak penalty as this reflects F I G U R E 3 Voronoi tesselation of the Trujillo trial, with each tile representing a house site. Also shown are the choropleth maps of the depth (left) and disc (right) measures of surroundedness. Intervention regions are those inside the thick black boundary the prior belief, in cluster randomized trials, that there is dependence within clusters and, for this particular trial, dependence within pairs. To obtain the posterior of the total intervention effect , as described in Section 3, we draw an independent sample of size 10 000 form the corresponding posterior marginals. All the orthogonalized models, as described in Section 2.4.4, were constructed by transforming the full ICAR design matrix derived from the neighborhood structure of the Voronoi tesselation. Figure 4 shows 95% posterior density intervals (quantile based) of the exponential of the total intervention effect (Tint), fitting the data with the standard model (2), the standard model (2) but with orthogonal spatial random effects and with the extended model (1) with the depth and disc measures of surroundedness (the latter with radius 200 m and both for the three levels of prior information). All the posterior intervals, apart from the standard model, cross the null value of 1, with medians on the beneficial side. Spatial dependence seems to be important. As pointed out in Section 2.4.1, orthogonalization only intends to increase the estimated variance of the fixed effects without major changes in the point estimates. In the extended model, most of the extra estimated variability is attributed to the added covariates (the surroundedness measure to model indirect effects), which are ignored in the standard models. Including indirect effects (ie, generalizing the standard model into the extended one) makes the intervals slightly wider when using the depth measure and much wider when using the disc measure of surroundedness. In other words, the difference in interval length implies, in this case, a strong indication of the presence of disc indirect effects but a weak indication of the presence of depth indirect effects.
The difference in interval length between extended depth and disc models is mainly due to the nature of the surroundedness measure itself. As shown in Figure 3, the range of values of the disc measure are comparable to those of the depth measure, but are spread more widely across the central part of the Trujillo area. In these terms, this area of the trial resembles the Crater configuration in F I G U R E 4 Note that 95% posterior density intervals for the exponential of the total intervention effect (Tint) in models fitted to data from the Trujillo trial. The points are the posterior medians. The upper three intervals are from models with the disc measure of surroundedness (radius 200 m), and the middle three for models with depth. Each measure was fitted with priors with strong, medium, or weak penalties. The bottom two intervals are from the standard model with and without spatial dependence Figure 1. As opposed to the models with depth measure, the effect of the prior penalty seems to be larger in the models with the disc measure, in which a stronger penalty implies a wider posterior interval for the total intervention effect.
Overall, in terms of posterior probability, it is clear that models with the depth measure provide stronger posterior evidence of a beneficial effect of the intervention as a whole. Figures 5 and 6 show 95% posterior intervals for the pairwise contralateral and ipsilateral effects, Sind 0 and Sind 1 , respectively, and for the same depth and disc models. As mentioned above, there is strong posterior evidence, in particular from the disc model, of an ipsilateral effect, that is, of locations in the intervention arm benefiting from being surrounded by other intervention locations, which was not a conclusion of the original analysis.

DISCUSSION
Geographically demarcated clusters are the most common choice for cluster-randomized trials (Hayes and Moulton, 2009). At the design stage, the clusters should ideally be sufficiently far apart to render spatial correlation negligible, as required by the standard analytic assumption of between-cluster independence. However, this may well be impractical, because the scale of spatial correlation can be up to tens of kilometers (Alexander et al., 2003;Clemens et al., 2006) and, in fact, may well be unknown at the design stage. Spatial aspects have been largely overlooked in the analysis of CRTs, beyond some incomplete ad hoc approaches (Benjamin-Chung et al., 2017;Jarvis et al., 2017). We have distinguished two possible consequences of spatial proximity in a CRT. Spillover dependence refers F I G U R E 5 Posterior 95% intervals for the exponential of the pairwise contralateral effect, for different models fitted to data from the Trujillo trial. The upper three intervals are from models with the disc measure of surroundedness (radius 200 m), and the lower three for models with depth. Each measure was fitted with strong, medium, or weak prior distributions to between-cluster dependence in individual outcomes, and spillover indirect effects increase or decrease the intervention effect, depending on the degree of surroundedness of individuals by other individuals in the intervention arm. We propose a new spatial regression model that incorporates both these phenomena, and whose intervention effects have a marginal interpretation (as opposed to being cluster-specific), as well as being well defined even when the proposed model is misspecified. We estimate, and adjust for, spillover indirect effects via a novel application of the statistical concept of depth (Tukey, 1975), and also using a simple measure of surroundedness based on the number of intervention individuals within a certain radius. The framework assumes spatial point data, although both surroundedness measures proposed can also be defined for areal data. For example, provided the centroids are good representatives of each area, for example, each lies within the area to it refers, we can use them to compute either of the two surroundedness measures. Minor modifications would also be required (eg, defining balls instead of discs and planes instead of lines) to define surroundedness in three-dimensional space (ie, including elevation).
The methodology proposed can be easily adapted to deal with total intervention effects defined in a different manner. For example, we can define the total intervention effect as a weighted average of expected rates at the cluster level (Hayes and Moulton, 2009) or as a simple average of all the pairwise intervention effects. However, the latter case would give unduly importance to large pairwise intervention effects that is not usually the objective in cluster randomized trials. CRT researchers would define total effects in such a way that, even though there are a few large pairwise effects, the intervention may still be considered effective as long as the total effects are small. Intervention effects may also be defined as geometrical averages instead F I G U R E 6 Posterior 95% intervals for the exponential of the pairwise ipsilateral effect, for different models fitted to data from the Trujillo trial. Labeling is as for Figure 5 of traditional arithmetic averages (see Hayes and Moulton, 2009, p.152).
In the Trujillo trial, there is some evidence that locations benefited indirectly from surrounding locations receiving the intervention. Despite being randomized, the clusters in each arm are distributed in just two contiguous regions. Restricted randomization could be used to achieve an allocation that addresses this or other concerns (see McCann et al. (2018) for one such proposed method). An alternative criterion could relate to geographical balance, in the intuitive sense of avoiding large contiguous regions of clusters from a single arm, although this notion may be hard to quantify. For example, despite lying in few contiguous regions, 17 of the 18 Trujillo clusters border at least one in the opposite arm. One approach may be to impose a minimum proportion of between-cluster boundaries that separate clusters in different arms, hence tending to reduce the size of contiguous single-arm regions.
We hope these methods will prove useful in those trials in which spatial correlation is inherent. Such analyses would provide more robust estimates of intervention effects, as well as quantifying the contribution of indirect effects, which are likely to be of interest when deciding whether to implement interventions at full scale.

A C K N O W L E D G M E N T S
This work was supported by the United Kingdom Medical Research Council (MRC) via Grant Reference G7508177, also by MRC and the UK Department for International Development (DFID), MRC Grant References MR/K012126/1 and MR/R010161/1, awards jointly funded by the MRC and DFID under the MRC/DFID Concordat agreement and also part of the EDCTP2 programme supported by the European Union. We also thank Christian Bottomley for helpful suggestions, Valeria Simoncini for guidance on iterative projection algorithms, and anonymous referees for useful suggestions.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings in this paper are available from the corresponding author upon reasonable request.