clusters. The outcome can be modeled as a function of exposure such that Y$_{i}$ = µ$_{i}$ ( β , A$_{i}$ ) = h ( β$_{0}$ + β$_{A}$A$_{i}$ ) , with h a link function. The general form of a classical GEE is U ( β ) = ∑ m i $_{=1}$D T $_{i}$V − 1 i { Y$_{i}$ − µ$_{i}$ ( β , A$_{i}$ ) } , where D$_{i}$ = ∂ µ$_{i}$ ( β ,A$_{i}$ ) ∂ β T is the design matrix, V$_{i}$ is the covariance matrix equal to φ R 1 / 2 i C ( α ) R 1 / 2 i , R$_{i}$ a diagonal matrix with elements var( Y$_{ij}$ ) , φ is the dispersion parameter, and C ( α ) is the “working” correlation structure with non-diagonal terms α . Parameters are estimated by setting U ( β ) to 0. Because our goal is to estimate the effect of exposure, our causal parameter of interest β$_{A}$ is the difference in average probability of being affected by the contagion process between the exposed and unexposed clusters: $_{E}$( Y | A = 1) −$_{E}$ ( Y | A = 0) . We assume each individual to have a set of P cluster-level and individual-level covariates X$_{ij}$ = ( X (1) ij , ..., X ( P ) ij ) $^{⊤}$, with all covariates represented compactly as X$_{i}$ = ( X$_{i}$$_{1}$ , ..., X$_{in}$$_{i}$ ) $^{⊤}$. Some of the covariates may relate to the network in which the individual is embedded and are described in the previous subsection. Accounting for this information is important in the estimation process. To account for non-randomized exposure A , it is possible to account for a different probability of exposure between clusters using a propensity score (PS) and to use classical inverse-probability weighting approaches for estimation of exposure effects. The propensity score is an arbitrary function g$_{i}$ ( X$_{i}$ , η$_{G}$ ) . The estimation procedure is unbiased when it exactly models P ( A$_{i}$ = a | X$_{i}$ ) . It is often fit using a logistic regression model if the cluster exposure status is binary and η$_{G}$ are nuisance parameters resulting from this estimation. To reduce confounding bias and improve efficiency of estimation, it is possible to introduce covariate adjustment to the standard GEE framework by augmenting the GEE itself. This requires specification of an outcome model (OM) B$_{i}$ ( X$_{i}$ , A$_{i}$ = a, η$_{B}$ ) = [ B$_{ij}$ ( X$_{i}$ , A$_{i}$ = a, η$_{B}$ )]$_{j}$$_{=1}$$_{,...,n}$$_{i}$ , which is an arbitrary function of X$_{i}$ given for each exposure level, where η$_{B}$ are nuisance parameters to be estimated. The estimation is unbiased and most efficient if the OM models the probability of the outcome of interest given baseline covariates$_{E}$ ( Y$_{ij}$ | X$_{i}$ , A$_{i}$ = a ) . The estimation is said doubly robust as either the OM or the PS (but not both) need to be set to a specific quantity to ensure unbiased estimation. In this case, we will call the OM or the PS “correctly specified” as it corresponds to the true data generation process. Equation 1 incorporates these additional terms in the estimating function. The n$_{i}$ × n$_{i}$ matrix G$_{i}$ ( X$_{i}$ , A$_{i}$, η$_{G}$ ) = diag [ A$_{i}$ g$_{i}$ ( X$_{i}$ ,η$_{G}$ ) + 1 − A$_{i}$ 1 − g$_{i}$ ( X$_{i}$ ,η$_{G}$ ) ] j =1 ,...,n$_{i}$ is an inverse probability weighting matrix specific to each cluster. The variance of β is estimated by using an empirical or “sandwich” estimator, which is also robust in the sense that it provides valid standard errors even when the assumed covariance structure is not correct. Given the OM and PS, the exposure effect is represented by the vector of coefficients β . 5 0 = m ∑ i =1 [ D T $_{i}$V − 1 i G$_{i}$ ( X$_{i}$ , A$_{i}$, η$_{G}$ ) ( Y$_{i}$ − B$_{i}$ ( X$_{i}$ , A$_{i}$, η$_{B}$ )) − ∑ a =0 , 1 D T i ( a ) V − 1 i ( a ) $^{(}$B$_{i}$ ( X$_{i}$ , A$_{i}$ = a, η$_{B}$ ) − µ$_{i}$ ( β , A$_{i}$ = a ) ] (1)