<figure><img src="image_1.png" /></figure>

HAL
open science

# Estimating an adjusted risk difference in a cluster randomized trial with individual-level analyses

Jules Antoine Pereira Macedo, Bruno Giraudeau

► To cite this version:

Jules Antoine Pereira Macedo, Bruno Giraudeau. Estimating an adjusted risk difference in a cluster randomized trial with individual-level analyses. Statistical Methods in Medical Research, 2025, 34 (1), pp.26-44. 10.1177/09622802241293783. hal-04952188

HAL Id: hal-04952188

https://hal.science/hal-04952188v1

Submitted on 17 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Estimating an adjusted risk difference in a cluster randomized trial with individual-level analyses

Jules Antoine Pereira Macedo¹ Bruno Giraudeau¹,², for the ESCIENT collaborators†

1. Université de Tours, Université de Nantes, INSERM, SPHERE U1246, Tours, France

2. INSERM CIC1415, CHRU de Tours, Tours, France

†. See “Acknowledgements” section.

Correspondence: Jules Antoine Pereira Macedo, Université de Tours, INSERM, SPHERE U1246, Tours, France.

Email: [pereiramacedo@univ-tours.fr](mailto:pereiramacedo@univ-tours.fr)

## Abstract

In cluster randomized trials (CRTs) with a binary outcome, intervention effects are usually reported as odds ratios, but the CONSORT statement advocates reporting both a relative and an absolute intervention effect. With a simulation study, we assessed several methods to estimate a risk difference (RD) in the framework of a CRT with adjustment on both individual- and cluster-level covariates. We considered both a conditional approach (with the generalized linear mixed model [GLMM]) and a marginal approach (with the generalized estimating equation [GEE]). For both approaches, we considered the Gaussian, binomial and Poisson distributions. When considering the binomial or Poisson distribution, we used the g-computation method to estimate the RD. Convergence problems were observed with the GEE approach, especially with low intra-cluster coefficient correlation values, small number of clusters, small mean cluster size, high number of covariates and prevalences close to 0. All methods reported no bias. The Gaussian distribution with both approaches and binomial and Poisson distributions with the GEE approach had satisfactory results in estimating the standard error. Results for type I error and coverage rates were better with the GEE than GLMM approach. We recommend using the Gaussian distribution because of its ease of use (the RD is estimated in one step only). The GEE approach should be preferred and replaced with the GLMM approach in cases of convergence problems.

Keywords: binary outcome, risk difference, cluster randomized trial, generalized estimating equations, generalized linear mixed models, g-computation
# 1. Introduction

In a cluster randomized trial (CRT), randomization units are groups of individuals, called clusters, rather than individuals themselves. When the outcome is binary, the CONSOLIDATED STANDARDS OF REPORTING TRIALS (CONSORT)¹,² statement recommends reporting both relative and absolute treatment effects by reporting both an odds ratio (OR) and/or a relative risk (RR) value as well as a risk difference (RD). Indeed, relative treatment effects often provide an overoptimistic view of the treatment effect.³⁻⁵ However, a recent systematic review of CRTs with binary primary outcomes³ showed a lack of compliance with these recommendations: in only 18% (13/73) of the trials were reported both a relative and an absolute treatment effect while in 63% (43/73) of the trials was reported only a relative treatment effect. Indeed, binary outcomes are classically analyzed with a logistic model, which allows for expressing intervention effects as ORs rather than as RDs.

Nevertheless, there are methods for estimating an RD, and simulation studies⁶⁻⁸ have been conducted. Ukoumunne et al.⁶ assessed the performance of a generalized estimating equation (GEE) approach using an identity link function to estimate an RD in a CRT. However, the authors did not consider any form of adjustment on covariates although CRTs are known to be prone to recruitment bias.⁹ Gehrmann et al.⁷ proposed estimating an adjusted RD considering one covariate with and without confounding, but with an individually randomized trial. Pedroza and Truong⁸ compared 6 methods using the GEE approach with different distributions and link functions, but they limited their study to multicenter individually randomized trials and considered only one balanced individual-level baseline covariate.

To our knowledge, no study has considered both individual- and cluster-level covariates in the framework of a CRT, although both types of variables may exist, with possible confounding for individual-level covariates. In this study, we investigated such situations and compared 6 methods to estimate an RD. We considered both generalized linear mixed model (GLMM) and GEE approaches, and for each approach, 3 probability distributions: Gaussian, binomial and Poisson.

# 2. Methods

We report our methods section following the guidelines, to generate data and reporting results of a simulation study, using the Aims, Data-generating mechanisms, Estimands, Methods, and Performances measures (ADEMP) structure as recommended by Morris et al.¹⁰

We considered a 2 parallel-arm CRT with a binary outcome. Individual and cluster-level covariates were generated, with individual-level covariates defined as confounding factors.

## 2.1. Simulation plan

### 2.1.1. Number of clusters and cluster size

Let *i* be the treatment arm index, with *i* = 0 for the control arm and *i* = 1 for the intervention arm. We considered a balanced design with randomization of 2*K* clusters in a 1:1 ratio with *j* being the cluster index within arm *i*, with *j* ∈ {1, ..., *k*}.
We denoted $n_{ij}$, the cluster size of the $j^{th}$ cluster in the $i^{th}$ arm, and $N$, the total number of individuals as $N = \sum_{i=0}^{1} \sum_{j=1}^{k} n_{ij}$. Cluster sizes were generated from a negative binomial distribution with parameters $n$, the mean cluster size, considering different mean cluster sizes when we had small and high number of clusters,¹¹ and $v$, the variance of the cluster size. Also, we took into account cluster size variability expressed by the coefficient of variation, $cv$, of the cluster size.¹² Thus, both $n$ and $cv$ were specified and $v$ was deduced as $v = cv^2n^2$. We set $cv$ at 0.6.¹²

### 2.1.2. Outcome

We denoted $Y_{ijl}$ as the outcome associated with the $l^{th}$ subject in the $j^{th}$ cluster of the $i^{th}$ arm. Outcomes were generated according to the Lunn and Davis¹³ approach as:

$$Y_{ijl} = (1 - U_{ijl})V_{ijl} + U_{ijl}Z_{ij} \quad (1)$$

where $U_{ijl}$, $V_{ijl}$, $Z_{ij}$ are independent Bernoulli variables defined as

$U_{ijl} \sim \text{Bern}(\sqrt{\rho})$, $V_{ijl} \sim \text{Bern}(p_i)$, $Z_{ij} \sim \text{Bern}(p_i)$ where $\rho$ is the intra-cluster correlation coefficient (ICC) and $p_i$ the prevalence in arm $i$. $Y_{ijl}$ is then Bernoulli distributed such that $Y_{ijl} \sim \text{Bern}(p_i)$, and we moreover have:

$$corr(Y_{ijl_1}, Y_{ijl_2}) = \rho, \quad (l_1 \neq l_2)$$

### 2.1.3. Covariates

We simulated covariates at both individual- and cluster-levels. The association between the outcome and covariates was expressed with a priori-specified ORs. Covariates were generated using the Jiang et al. approach (algorithm 4).¹⁴

#### 2.1.3.1. Individual-level covariates

Let $X_{ijl}^{(m)}$ be the $m^{th}$ individual-level covariate of the $l^{th}$ individual of the $j^{th}$ cluster of the $i^{th}$ arm. $X^{(m)}$ are binary covariates defined as:

$$X_{ijl}^{(m)} = A_{ijl}^{(m)} W_{ijl}^{(m)}$$

with

- $A_{ijl}^{(m)} \sim \text{Bern}(\alpha_m)$, with $\alpha_m = \frac{q_m}{p_i}$ where $q_m$ is the prevalence of the $m^{th}$ covariate

- $W_{ijl}^{(m)} = (1 - R_{ijl}^{(m)})G_{ijl}^{(m)} + R_{ijl}^{(m)}Y_{ijl}$, with $R_{ijl}^{(m)} \sim \text{Bern}(r_{i,m})$ with $r_{i,m} = \rho_m \frac{\sqrt{(1-q_m)}}{\sqrt{\alpha_m(1-p_i)}}$ and $G_{ijl}^{(m)} \sim \text{Bern}(p_i)$,

- $A_{ijl}^{(m)}$ and $W_{ijl}^{(m)}$ are independent. $R_{ijl}^{(m)}, G_{ijl}^{(m)}$ and $Y_{ijl}$ are independent.
According to this approach, $X^{(m)}$ is Bernoulli-distributed such that $X_{ijl}^{(m)} \sim Bern(q_m)$ (Appendices 1 and 2), and we also have $corr(X_{ijl}^{(m)}, Y_{ijl}) = \rho_m$, (Appendix 3).

Such a correlation between the covariate $X^{(m)}$ and the outcome $Y_{ijl}$ can also be expressed as an OR. Indeed, the calculation of the OR between $X^{(m)}$ and $Y_{ijl}$ is as follows (Appendix 4):

$$OR_m = \frac{(q_m r_{i,m} + p_i q_m - p_i q_m r_{i,m})(1 - q_m + q_m r_{i,m} - p_i + p_i q_m - p_i q_m r_{i,m})}{(q_m - q_m r_{i,m} - p_i q_m + p_i q_m r_{i,m})(p_i - p_i q_m - q_m r_{i,m} + p_i q_m r_{i,m})} \quad (2)$$

From a practical point of view, we could first a priori specify an OR (which is more convenient to specify than the correlation between 2 binary variables) and deduce the value of $\rho_m$ by resolving equation (2), such that the solution belongs to the interval

$$\max \left\{ -\left( \frac{p_1 q_m}{(1-p_1)(1-q_m)} \right)^{\frac{1}{2}}, -\left( \frac{(1-p_1)(1-q_m)}{p_1 q_m} \right)^{\frac{1}{2}}, -\left( \frac{p_0 q_m}{(1-p_0)(1-q_m)} \right)^{\frac{1}{2}}, -\left( \frac{(1-p_0)(1-q_m)}{p_0 q_m} \right)^{\frac{1}{2}} \right\}; \\ \min \left\{ -\left( \frac{p_1(1-q_m)}{(1-p_1)q_m} \right)^{\frac{1}{2}}, -\left( \frac{(1-p_1)q_m}{p_1(1-q_m)} \right)^{\frac{1}{2}}, -\left( \frac{p_0(1-q_m)}{(1-p_0)q_m} \right)^{\frac{1}{2}}, -\left( \frac{(1-p_0)q_m}{p_0(1-q_m)} \right)^{\frac{1}{2}} \right\}$$

as required for a correlation coefficient between 2 binary variables.¹⁵

The intra-cluster correlation for covariate $X^{(m)}$ then equals $corr(X_{ijl_1}^{(m)}, X_{ijl_2}^{(m)}) = \rho_m^2 \rho$.

Such a definition of $X^{(m)}$ implies that $q_m < \min(p_0, p_1)$ thus allowing for $\alpha_m \in [0,1]$. However, by considering the complementary covariate $\tilde{X}^{(m)}$ defined as $\tilde{X}^{(m)} = 1 - X^{(m)}$, we may obtain a new covariate with prevalence $\tilde{q}_m = 1 - q_m$, whose correlation with the outcome equals $corr(\tilde{X}^{(m)}, Y_{ijl}) = -\rho_m$. Thus, depending on the desired prevalence of the covariate, one may use $X^{(m)}$ or $\tilde{X}^{(m)}$.

### 2.1.3.2. Confounding in individual-level covariates

Confounding was introduced by using the Leyrat et al.¹⁶ approach. We generated a selection variable $D_{ijl}$, specifying whether the $l^{th}$ individual of the $j^{th}$ cluster of the $i^{th}$arm was finally included in the trial, with $D_{ijl} \sim Bern(e_{ijl}^{(D)})$. Hence $e_{ijl}^{(D)}$ corresponds to the individual propensity score to be included in the trial, depending on individual-level covariates and a cluster effect. It is defined as:

$$e_{ijl}^{(D)} = \frac{1}{1 + e^{-L_{ijl}}}, \quad (3)$$

with

$$L_{ijl} = \sum_{m=1}^{S} \beta_{(m),i} (X_{ijl}^{(m)} - E[X^{(m)}]) + \zeta_{ij}, \quad (4)$$

with $E[X^{(m)}]$ the expected value of the $m^{th}$ individual-level covariate (corresponding to the prevalence $q_m$, as defined previously) and $\beta_{(m),i}$ the parameter quantifying the association between the $m^{th}$ individual-level covariate with the inclusion status in the $i^{th}$ arm. $\zeta_{ij}$ is a cluster effect associated with an inclusion rate $\pi_i$ for the
i'th arm, with $\zeta_{ij} = \ln(\frac{\tau_{ij}}{1-\tau_{ij}})$, with $\tau_{ij}$ the inclusion rate for the j'th cluster of the i'th arm, defined as $\tau_{ij} \sim Beta(\frac{\pi_i(1-\rho_D)}{\rho_D}, \frac{(1-\pi_i)(1-\rho_D)}{\rho_D})$, where $\rho_D$ is an inclusion ICC. Parameters a priori set were $\beta_{(1),1} = 2.5, \beta_{(2),1} = 1.5, \beta_{(3),1} = 1$, and $\beta_{(1),0} = 1, \beta_{(2),0} = 1, \beta_{(3),0} = 1$ thus corresponding to the inclusion status not being associated with individual-level covariates in the control arm but associated in the experimental arm. By setting different values for $\beta_{(1)}, \beta_{(2)}$ and $\beta_{(3)}$, depending on the intervention arm, the distribution of the covariates will differ between intervention arms and as explained in the section "Individual-level covariates". Because these covariates are also associated with the outcome, this means that they are associated with both the outcome and the randomization arm, thus defining confounding. We also specified $\pi_1 = \pi_0 = 0.7$, corresponding to a 70% rate of inclusion in each arm and $\rho_D = 0.3$ as in Leyrat et al.¹⁶. Because inclusion rates $\pi_i$ are < 1, one must increase the mean cluster size, $n$, by a factor $\frac{1}{0.7}$ to obtain the required $n$ value.

### 2.1.3.3. Cluster-level covariates

Let $X^{(m)}$ be the m'th cluster-level covariate defined as:

$$X_{ij}^{(m)} = A_{ij}^{(m)} W_{ij}^{(m)}$$

with

- $A_{ij}^{(m)} \sim Bern(\alpha'_m)$, with $\alpha'_m = \frac{q'_m}{p_i}$ where $q'_m$ is the m'th cluster covariate prevalence

- $W_{ij}^{(m)} = (1 - R_{ij}^{(m)})G_{ij}^{(m)} + R_{ij}^{(m)}Z_{ij}$, where $R_{ij}^{(m)} \sim Bern(r'_m)$ $G_{ij}^{(m)} \sim Bern(p_i)$ and $Z_{ij} \sim Bern(p_i)$ as defined in equation (1).

- $A_{ijl}^{(m)}$ and $W_{ijl}^{(m)}$ are independent. $R_{ij}^{(m)}, G_{ij}^{(m)}$ and $Z_{ij}$ are independent.

According to this approach, $X^{(m)}$ is Bernoulli-distributed such that $X_{ij}^{(m)} \sim Bern(q'_m)$. One may specify an OR associating $X^{(m)}$ with the outcome, from which we can then derive the $r'_m$ parameter by solving the following equation, with the condition of obtaining a solution belonging to the [0,1] interval (Appendix 5):

$$OR'_{m} = \frac{(p_{i}q'_{m} + q'_{m}r'_{m}\sqrt{\rho} - p_{i}q'_{m}r'_{m}\sqrt{\rho})(1 - q'_{m} + q'_{m}r'_{m}\sqrt{\rho} - p_{i} + p_{i}q'_{m} - p_{i}q'_{m}r'_{m}\sqrt{\rho})}{(q'_{m} - q'_{m}r'_{m}\sqrt{\rho} - p_{i}q'_{m} + p_{i}q'_{m}r'_{m}\sqrt{\rho})(p_{i} - p_{i}q'_{m} - q'_{m}r'_{m}\sqrt{\rho} + p_{i}q'_{m}r'_{m}\sqrt{\rho})} \quad (5)$$

### 2.1.4. Simulation plan parameters

Values of the simulation parameters are in Table 1, for a total of 1744 scenarios. The parameters were selected from previous simulation studies on CRTs such as the studies of Leyrat et al.¹¹ and Thompson et al.¹⁷, whose parameter values were based on CRT reviews.

[Insert Table 1]
## 2.2. Analysis methods

We considered 2 different approaches to analyze data: the cluster-specific approach by using the GLMM¹⁸ and the population averaged approach by using the GEE.¹⁹ We then assessed 6 methods, with 3 distributions for each approach, Gaussian, binomial and Poisson, and their respective link functions (i.e., identity, logit and log).

When considering the Gaussian distribution with an identity link function, we directly obtained RD estimates. However, when considering binomial or Poisson distributions, we actually estimated log ORs or log RRs and therefore used the g-computation method to estimate the RD.

### 2.2.1. GLMM

Let us consider the following GLMMs:

$$g(\pi_{ijl}) = \alpha + \theta T_{ij} + \beta X_{ijl} + \beta' X'_{ij} + \mu_{ij}, \quad (6)$$

with $\pi_{ijl}$ the probability of success of the $l^{th}$ individual of the $j^{th}$ cluster of the $i^{th}$arm and $g(\cdot)$ the link function (identity, logit, or log), determined by the probability density function considered (Gaussian, binomial, or Poisson, respectively). Fixed effects are $\alpha$, for the intercept, $\theta$, for the intervention effect, $\beta = (\beta_1, ..., \beta_S)$ a S vector associated with individual-level covariates and $\beta' = (\beta'_1, ..., \beta'_S)$, a S' vector associated with cluster-level covariates. $T_{ij}$ is a binary variable associated with the treatment arm. Finally, $\mu_{ij}$ are random effects with $\mu_{ij} \sim N(0, \sigma^2)$. Both fixed and random effects coefficients were estimated with the restricted maximum likelihood method when considering the Gaussian distribution and maximum likelihood method when considering the binomial and Poisson distributions.

### 2.2.2. GEE

The GEE approach modeled the probability of success as:

$$g(\pi_{ijl}) = \alpha + \theta T_{ij} + \beta X_{ijl} + \beta' X'_{ij} \quad (7)$$

with the parameters corresponding to those used for the GLMM approach. An exchangeable correlation structure was considered, and estimates were obtained using the quasi-likelihood approach, whatever the distribution.

### 2.2.3. G-computation

G-computation was proposed by Robins.²⁰ The aim is to predict the probability of success of each individual in both situations (i.e., whether the individual receives the intervention [experimental arm] or not [control arm]). Hence, we must consider a counterfactual event, defined as the event in the hypothetical situation in which the individual would be in the other arm than the arm to which they had actually been allocated. Thus, for an individual allocated to the experimental arm, the counterfactual event would be the event that could have been observed if the individual had been allocated to the control arm and vice versa. G-computation has 3 main steps:
1) We estimated the probability of success by estimating the parameters with GLMM and GEE approaches, as expressed in equations (6) and (7).

2) For each individual, we predicted the probability of success and the counterfactual event using the parameter estimates obtained in the first step. We denoted $\hat{\pi}_{ijl}$, the estimation of the success probability of the l$^{th}$ individual of the j$^{th}$ cluster of the i$^{th}$arm and $\hat{\pi}_{ijl}^{(c)}$, the estimation of the counterfactual event. For GLMM (as for the GEE approach), we can denote these estimations as:

$$ \hat{\pi}_{ijl} = g^{-1}(\hat{\alpha} + \hat{\theta}T_{ij} + \hat{\beta}X_{ijl} + \hat{\beta}'X'_{ij}) \quad (8) $$

$$ \hat{\pi}_{ijl}^{(c)} = g^{-1}(\hat{\alpha} + \hat{\theta}|T_{ij} - 1| + \hat{\beta}X_{ijl} + \hat{\beta}'X'_{ij}) \quad (9) $$

As we can see in equations (8) and (9), we did not consider the random effect for the GLMM approach (see below in section “Software”).

3) We estimated individual RDs, $\hat{RD}_{ijl}$, as the difference in these probabilities as:

$$ \hat{RD}_{ijl} = \begin{cases} \hat{\pi}_{ijl} - \hat{\pi}_{ijl}^{(c)}, & \text{if } i = 1 \\ \hat{\pi}_{ijl}^{(c)} - \hat{\pi}_{ijl}, & \text{if } i = 0 \end{cases} \quad (10) $$

Thus, we obtained an estimation of the RD as the average of these individual RDs:

$$ \hat{RD} = \frac{1}{N} \sum_{i=0}^{1} \sum_{j=1}^{k} \sum_{l=1}^{n_{ij}} \hat{RD}_{ijl} \quad (11) $$

### 2.2.4. Standard error (SE)

When using the Gaussian distribution, the SE of the RD was directly estimated. When considering the binomial and Poisson distributions, we used the delta method²¹ to estimate SE. However, in CRTs, SEs are not well estimated, specifically in CRTs with a small number of clusters. Hemming and Taljaard²² advocated using a small sample size correction and the t-distribution to correct the variance under-estimation due to the CRT design. Li and Redden²³,²⁴ tested multiple small sample size corrections for GLMM and GEE approaches. For GLMM, this corresponds to a correction of the denominator degree of freedom.²⁵⁻²⁷ According to their work and our simulation parameters, we chose the between-within²⁵ (B-W) method. Using a t-distribution with the B-W correction method, the degree of freedom then equals the number of clusters minus the number of fixed-effect coefficients to be estimated [i.e., 2k - (S + S' + 2)]. For GEE, the variance-covariance of the regression coefficient estimators can be estimated by a sandwich estimator, using a correlation type matrix. This estimator is biased and under-estimates the variances and covariances. Some methods²⁸⁻³¹ have been proposed to improve the sandwich estimator; we chose the Fay and Graubard correction²⁸ and used a t-distribution with a degree of freedom equal to the number of clusters minus the number of parameters estimated.

### 2.2.5. Performance measures
We simulated $n_{sim} = 1000$ datasets for each of the 1744 scenarios and evaluated 5 performance measures. We chose 1000 iterations to allow for estimating a reasonable Monte Carlo SE of 0.689 for type I error rate and coverage rate performance measures.¹⁰ The 5 performance measures used in this simulation study are the convergence rate, bias, type I error rate, coverage rate, and finally the relative error of the SE rate. For the convergence rate, for each scenario and method, we determined the rate from the 1000 iterations for which the model was fitted. For bias, we estimated the difference between the a priori-specified RD and the estimator. For type I error rate (i.e., when considering scenarios with no treatment effect), we estimated the rate of iterations, among the convergent iterations, rejecting the null hypothesis of no treatment effect by determining whether the true value (RD = 0) was out of the estimated confidence interval. For coverage rate (i.e., considering only scenarios with a non-null treatment effect), we estimated the rate of iterations, among the convergent iterations, with the a priori-specified RD in the estimated confidence intervals. Finally, for relative error of the SE rate, we calculated the ratio between the average estimated SEs and the estimated empiric SE. Formulas to determine these performance measures are in Table 2. The type I error $\alpha$ was fixed at 5%.

[Insert Table 2]

### 2.2.6. Software

To perform our simulation study, we used R v4.1.1.³² We used the function *lmer* and *glmer* from the *lme4* package for GLMMs and compared the *glmer* function with the *glmmPQL* function from the *MASS* package. For predicting the individual success probabilities, contrary to the *glmmPQL* function, *glmer* did not consider the random effects using the *predict* function from the *stats* package. We did observe negligible differences in the results depending on whether the *glmer* or the *glmmPQL* function was used with the *predict* function (only the 5th decimal changed). Reported results are those obtained with the *glmer* function. We used the *gee* function from the *gee* package for the GEE approach and the *GEE.var.fg* function, which we modified so that it could account for variable cluster sizes, from the *geesmv* package for the small sample size correction of Fay and Graubard. We used the *calcl centre from the University of Orléans-Tours to parallelize the simulation and used the *parallel* package. We set 2 random seeds depending on scenarios with no treatment effect and scenarios with a non-null treatment effect. For the 972 scenarios with no treatment effect, we set *set.seed(37250)*. For the 972 scenarios with a non-null treatment effect, we set *set.seed(1995)*.

The code for the simulation study is available at

https://github.com/JulesPereiraMacedo/RD_estimation_CRT_indiv_analyses.

## 3. Results

Figures 1, 3, 4 and 6 report the overall performance of the methods at each level of the simulation parameters.

Supplementary figures S1 to S72 complete the information for specific parameter combinations. Figures 2 and 5 report the bias of all scenarios, with the combination of 3 parameters (number of clusters per arm, ICC and total number of covariates). Any specific combination of the simulation parameters can be found in the supplementary figures.
## 3.1. No treatment effect

### 3.1.1. Convergence rates

With the GLMM approach, we had good convergence rates whatever the distribution. Convergence failures were indeed anecdotal: less than 0.1% for the Gaussian distribution and less than 0.4% for the binomial and Poisson distributions (Supplementary Table 1). The GEE approach produced convergence problems with all 3 distributions. These problems increased when the ICC value decreased, the number of clusters per arm decreased, the number of covariates increased, the mean cluster size decreased, or the common prevalence decreased. We did not observe systematic differences according to distribution (Supplementary Table 1). Referring results are in supplementary file Table 1.

### 3.1.2. Type I error

We considered the interval $CI_{TIE} = [3%; 7%]$, which corresponds to a 99.6% interval for the distribution of the estimated type I error rate, considering that the Monte Carlo SE of the estimated type I error was 0.689 because the number of iterations specified was 1000. Figure 1 reports the proportion of scenarios for which the estimated type I error (TIE) rate was within or outside of the $CI_{TIE}$ interval. The methods were considered conservative if the type I error was under the 3% limit and liberal if over the 7% limit. Results were poor when using the GLMM approach, with a large proportion of scenarios for which this approach was conservative, whatever the probability distribution. Results were even worse with the Poisson distribution versus the 2 other distributions, with almost 1 in 2 scenarios for which the type I error rate was considered acceptable. Results with the GEE approach were much better, with no major differences by distribution, although the method based on the Gaussian distribution was slightly more conservative than with the other two distributions. This conservatism seemed to increase with both the GLMM and GEE approaches as the number of clusters per arm decreased, slightly when the number of covariates increased, and also when the mean cluster size decreased. The impact of the ICC differed between the two approaches: for the GLMM method, when the ICC decreased, the conservatism increased, whereas for the GEE method, it slightly increased when the ICC increased. Conservatism was also increased as the common prevalence decreased except for the GLMM approach when using the Poisson distribution, for which it increased. The GLMM and GEE approaches were more liberal as the number of clusters per arm decreased, the ICC increased, the number of covariates increased, the common prevalence increased, and the mean cluster size decreased. The higher the number of covariates, the lower the proportion of scenarios for which the type I error rate was between 3% and 7%. Actually, both the number of clusters and number of covariates affect the degree of freedom of the Student distribution. Thus, with a small number of clusters per arm, such as $k = 5$, the degrees of freedom can be 3 (if 3 individual-level covariates and 2 cluster-level ones), which may explain the observed results.

[Insert Figure1]

### 3.1.3. Bias
In Figures 2, S13 to S24, the bias estimated increased when the ICC increased, especially with an ICC value of 0.2, and the number of clusters per arm decreased. The maximum absolute value of bias was 0.012, which implies that all methods were free of bias. These results were expected because all our analyses were adjusted on the confounding variables.

[Insert Figure 2]

### 3.1.4. Relative bias of the SE

We looked at the relative bias of the SE with the limit values -10% and +10%: the SE was considered under-estimated if under the -10% limit and over-estimated if over the +10% limit. For the 2 methods based on the Gaussian distribution and the methods using the GEE approach (with a binomial or Poisson distribution) the relative bias of the SE was acceptable in most scenarios (Figure 3). With these 4 methods, there was some form of over-estimation of the SE notably for low ICC values, small number of clusters per arm, small mean cluster size, and high common prevalences. The number of covariates had no impact. Of note, under-estimation of the SE was anecdotal, thus illustrating the efficacy of the Fay and Graubard correction methods for a small number of clusters when using a GEE approach. For the 2 last methods (GLMM with the binomial or Poisson distribution), the SE was under-estimated in a large proportion of scenarios. This was amplified as the ICC increased, the number of clusters decreased and, to some extent, the number of covariates increased. In some situations, the SE was under-estimated in half of the scenarios.

[Insert Figure 3]

## 3.2. Non-null treatment effect

### 3.2.1. Convergence rates

Again, the GLMM approach had good convergence rates whatever the distribution. The failure rates of the GLMM convergence were slightly higher than in scenarios with no treatment effect, but they still remained anecdotal: no convergence failure with the Gaussian distribution, a maximal value of 1.5% with the binomial distribution and 1.2% with the Poisson distribution (Supplementary Table 2). With the GEE approach, convergence problems occurred with all 3 distributions. These problems of convergence were severe: at least 20% to 25% in some scenarios. They increased when the ICC value decreased, the number of clusters per arm decreased, the number of covariates increased, the mean cluster size decreased and the prevalences in the 2 groups decreased. We found no systematic differences by distribution. These results were similar to those observed in the convergence rates analysis for scenarios with no treatment effect. These convergence problems were already observed in CRT simulation studies using the GEE approach with an exchangeable correlation matrix, especially when the number of clusters was small<sup>6,17</sup> and with a small degree of clustering<sup>33,34</sup>. Referring results are in supplementary file Table 2.

### 3.2.2. Coverage rates
We considered the interval $CI_{CR} = [93\%; 97\%]$, which corresponds to a 99.6% interval for the distribution of the estimated coverage rate, considering that the Monte Carlo SE of the estimated coverage rate was 0.689 because the number of iterations specified was 1000. Figure 4 reports the proportion of scenarios for which the estimated coverage rate (CR) was in the $CI_{CR}$ interval. The methods were considered conservative if the CR was over the 97% limit and liberal if under the 93% limit. With the GLMM approach, for a low proportion of scenarios, the coverage rate could be considered acceptable, whatever the distribution used. Of note, the method based on the Gaussian distribution was consistently conservative, except for scenarios with an ICC equal to 0.2, for which it was liberal, whereas the methods based on the binomial and Poisson distributions were liberal for some scenarios and conservative for others. The results from the GEE approach were much better, with no major difference by distribution. The conservatism seemed to increase when the number of clusters per arm decreased, the ICC decreased, the number of covariates decreased, the prevalence increased and, the mean cluster size decreased. We observed this conservatism for the methods using the GLMM approach. For the GEE approach, we observed a low proportion of scenarios for which the methods were conservative, when the number of clusters per arm was low. Liberalism increased when the number of clusters per arm decreased, the ICC increased, when there were covariates, and in situations of pairs of prevalence of 0.5/0.4 and 0.1/0.05, and when the mean cluster size increased. These observations share similarities with those in the “Type I error” section, as we could expect. The similarity in results obtained with the GEE approach, whatever the distribution considered, were also observed in previous studies.⁶⁻⁸ Otherwise, the conservatism we observed if the number of clusters per arm was small agrees with the Ukoumune et al. results from a simulation study conducted in the framework of a CRT.⁶ Finally, as already mentioned in the “Type I error” section, the degrees of freedom of the Student distribution used can explain the conservatism of our methods in some scenarios.

[Insert Figure 4]

### 3.2.3. Bias

In Figure 5, and for all scenarios (S49 to S60), we observed no bias with any methods, except for the GEE approach using the Poisson distribution when the ICC increased to 0.2. Our results agree with those observed in trials assessing a non-null intervention effect and considering clustering,⁶,⁸,¹⁷ even though these latter studies did not consider confounding. Of note, these observed results were expected because all our analyses were adjusted on the confounding variables.

[Insert Figure 5]

### 3.2.4. Relative bias of the SE

Again, we considered the -10% to +10% interval to identify scenarios for which the SE estimation could be considered acceptable. As for scenarios with no treatment effect, for scenarios with a no-null treatment effect, the 2 methods based on the Gaussian distribution and the 2 other methods that used the GEE approach (with a binomial or Poisson distribution) produced a high proportion of scenarios for which the relative bias of the SE was acceptable (Figure 6). Conversely, when the GLMM approach was used with the binomial or Poisson distribution, the SE was poorly estimated: the under-estimation was more pronounced than for all 4 other
methods in all scenarios and was more often under-estimated when the binomial distribution was considered.

Also, the over-estimation was more pronounced in all scenarios when the Poisson distribution was considered.

The influence of the ICC was modest, except when the ICC was very low ($\rho$ = 0.001), which probably illustrates a difficulty in correctly estimating the SE when there is nearly no correlation in data. Low number of clusters and small mean cluster size were also associated to difficulties in correctly estimating the SE. With the GEE approach, results were even more pronounced when using the Poisson distribution versus the 2 other distributions, perhaps because of poor performance of the Fay and Graubard correction for a small number of clusters with the Poisson distribution. The results of this section are consistent with the results of the “Relative bias of the SE” section obtained for scenarios with no treatment effect.

[Insert Figure 6]

## 4. Motivating example

The Ambroisie trial³⁵ was a CRT aimed at demonstrating that continued enteral nutrition until extubation in patients in an intensive care unit (ICU) was not inferior to fasting before extubation. Twenty-two ICUs in France were randomized 1:1 to the 2 treatments. The primary outcome was the rate of extubation failure defined as a composite criterion: reintubation or death within the 7 days after extubation. A total of 1130 patients were included and, in the intention-to-treat population, the RD difference was estimated at -0.4% (95% CI -5.2% to 4.5%); the non-inferiority margin was a priori specified as 10%. The RD was estimated by the difference in proportions between both groups. No adjustments were performed. The CI was estimated by correcting the variance using a clustering correction factor to consider the inflation in variance of each proportion in each arm due to clustering.³⁶

We selected individual- and cluster-level covariates with possible individual-level confounding covariates. For individual-level covariates, we considered the obesity status, the use of corticosteroids (both when used to prevent the “laryngeal oedema” before extubation and also 24 hours before extubation for any other indication), whether the patient had a chronic respiratory disease or not (chronic obstructive pulmonary disease, asthma, obesity-hypoventilation syndrome and restrictive respiratory disease), and, finally whether the patient had heart disease or not (left ventricular dysfunction, chronic atrial fibrillation, history of pulmonary edema, documented ischemic cardiomyopathy and/or severe valvular disease). We also considered the university status of the hospital as a cluster-level covariate. We adjusted for all these variables. Because these covariates had missing data (11.5% of the baseline data), we used multiple imputation with the predictive mean matching method, which was performed 10 times on the dataset with the covariates mentioned above before the binarization. Results are in Figure 7. We observed similar point estimations of the RD as the RD reported in the study. The 95% CIs were slightly wider with the GEE than GLMM approach. Whatever the method used, non-inferiority would have been concluded.

[Insert Figure 7]

## 5. Discussion
We compared 6 methods to estimate an RD in a CRT with adjustment on both cluster-level covariates and individual-level confounding covariates. We observed consistent results, whether scenarios had an intervention effect or not. We did not observe any form of bias, whatever the method used. The relative bias of the SE was acceptable in a large proportion of scenarios with the GLMM and GEE approaches when using Gaussian distributions and also for the 2 other GEE methods (binomial and Poisson distributions), which shows the efficiency of the Fay and Graubard correction for small number of clusters using the GEE approach.

Nevertheless, although the type I error and coverage rates were satisfactory when using the GEE approach whatever the distribution, they were not acceptable in many scenarios when using the GLMM approach even with the Gaussian distribution. One explanation may be that the between-within correction method used for small number of clusters was not the most appropriate. Indeed, in some scenarios with a small number of clusters and a high number of covariates, the corrected degrees of freedom may be as small as 3, which may lead to 2 wider CIs and thus explain the conservatism of GLMM approach. We observed better results for performance measures when the number of clusters per arm increased. For type I error and coverage rates, results were satisfactory when the number of clusters per arm was 30 or more. Hence, the GEE approach offers better properties than the GLMM approach, although we must keep in mind the convergence problems, which were not observed with the GLMM approach. As complementary analyses, for scenarios with < 80% convergence, we looked at the performance of the composite method in which the GEE approach is used first and replaced by the GLMM approach in case there is non-convergence. Both approaches were used with the Gaussian distribution.

We observed (data not shown) that the composite method allowed for recovering the performance of the GEE approach for both type I error and coverage rates and over-performed the GLMM approach. As for the relative error of the SE rate, the composite method overperformed both the GEE and GLMM approaches.

A limitation of this work is that we addressed only the participant-average estimand, with the marginal and conditional approach, as explained by Kahan et al.³⁷ In doing so, in each method, all participants are equally weighted, but we could have considered cluster-average estimands for which the clusters are equally weighted. We indeed did not compare the participant-average and cluster-average estimand. Another limitation is that we used only one correction method for small number of clusters for each approach despite other available methods. Finally, we did not consider the random effect neither in the g-computation nor when using the delta method for the GLMM approach.

Beyond these limitations, a strength of this simulation study is that we considered situations with more than one covariate, as is often done. To our knowledge, no simulation study in the framework of a CRT considered both individual- and cluster-level covariates. We also considered confounding individual-level covariates, which may occur in CRTs and is seldom considered.

## 6. Conclusion

We recommend use of a method based on the Gaussian distribution with the GEE approach because of its ease of use (an RD can be estimated in a single step, directly from a model), rather than the binomial or Poisson distribution, which requires using both the g-computation and delta methods. Moreover, we advocate using the GEE approach with an exchangeable correlation matrix because of its performance in estimating the SE, which translates to good type I error and coverage rates. Finally, for cases of convergence problems with the GEE
approach, we suggest fitting a GLMM with a Gaussian distribution. Also, these convergence problems and poor
performance measures for scenarios with 10 or fewer clusters per arm highlighted that better approaches such as
Targeted Maximum Likelihood Estimation³⁸,³⁹, are needed in the CRT framework to estimate RD.

7. Acknowledgements

†The authors thank Nelly Agrinier, Laetitia Minary, Joëlle Kivits and Marie Buzzi working on the ESCIENT project. The authors thank the scientific board: Monica Taljaard, Elizabeth Turner, Sandra Eldridge and Agnès Caille.

The authors thank Mickaël Landais (Le Mans Hospital) and Stephan Ehrmann (Tours University Hospital) for sharing the data of the AMBROISIE trial. The authors also thank Laura Smales for English proofreading. The authors benefitted from the use of the cluster CaSciModOT at the Datacenter Centre-Val de Loire.

8. Declaration of conflicting interests

The authors declares that there is no conflicts of interest.

9. Funding

This work was supported by the French National Research Agency [ANR-21-CE36-0007-01-ESCENT].
References

1. Moher D, Hopewell S, Schulz KF, et al. CONSORT 2010 Explanation and Elaboration: updated guidelines for reporting parallel group randomised trials. *BMJ* 2010; 340: c869.

2. Campbell MK, Piaggio G, Elbourne DR, et al. Consort 2010 statement: extension to cluster randomised trials. *BMJ* 2012; 345: e5661.

3. Turner EL, Platt AC, Gallis JA, et al. Completeness of reporting and risks of overstating impact in cluster randomised trials: a systematic review. *Lancet Glob Health* 2021; 9: e1163–e1168.

4. Malenka DJ, Baron JA, Johansen S, et al. The framing effect of relative and absolute risk. *J Gen Intern Med* 1993; 8: 543–548.

5. Gallis JA, Turner EL. Relative Measures of Association for Binary Outcomes: Challenges and Recommendations for the Global Health Researcher. 2019; 85: 137.

6. Ukoumunne OC, Forbes AB, Carlin JB, et al. Comparison of the risk difference, risk ratio and odds ratio scales for quantifying the unadjusted intervention effect in cluster randomized trials. *Stat Med* 2008; 27: 5143–5155.

7. Gehrmann U, Kuss O, Wellmann J, et al. Logistic regression was preferred to estimate risk differences and numbers needed to be exposed adjusted for covariates. *J Clin Epidemiol* 2010; 63: 1223–1231.

8. Pedroza C, Truong VT. Performance of models for estimating absolute risk difference in multicenter trials with binary outcome. *BMC Medical Res Methodol* 2016; 16: 113.

9. Higgins JP, Eldridge S, Li T. Chapter 23: Including variants on randomized trials. In: *Cochrane Handbook for Systematic Reviews of Interventions*. 2023.

10. Morris TP, White IR, Crowther MJ. Using simulation studies to evaluate statistical methods. *Stat Med* 2019; 38: 2074–2102.

11. Leyrat C, Morgan KE, Leurent B, et al. Cluster randomized trials with a small number of clusters: which analyses should be used? *Int J Epidemiol* 2018; 47: 321–331.

12. Eldridge SM, Ashby D, Kerry S. Sample size for cluster randomized trials: effect of coefficient of variation of cluster size and analysis method. *Int J Epidemiol* 2006; 35: 1292–1300.

13. Lunn AD, Davies SJ. A Note on Generating Correlated Binary Variables. *Biometrika* 1998; 85: 487–490.

14. Jiang W, Song S, Hou L, et al. A Set of Efficient Methods to Generate High-Dimensional Binary Data With Specified Correlation Structures. *Am Stat* 2021; 75: 310–322.

15. Prentice RL. Correlated Binary Regression with Covariates Specific to Each Binary Observation. *Biometrics* 1988; 44: 1033–1048.

16. Leyrat C, Caille A, Donner A, et al. Propensity scores used for analysis of cluster randomized trials with selection bias: a simulation study. *Stat Med* 2013; 32: 3357–3372.

17. Thompson JA, Leyrat C, Fielding KL, et al. Cluster randomised trials with a binary outcome and a small number of clusters: comparison of individual and cluster level analysis method. *BMC Med Res Methodol* 2022; 22: 222.

18. Breslow NE, Clayton DG. Approximate Inference in Generalized Linear Mixed Models. *J Am Stat Assoc* 1993; 88: 9–25.

19. Zeger SL, Liang K-Y, Albert PS. Models for Longitudinal Data: A Generalized Estimating Equation Approach. *Biometrics* 1988; 44: 1049–1060.
20. Robins J. A new approach to causal inference in mortality studies with a sustained exposure period—
application to control of the healthy worker survivor effect. *Math Model* 1986; 7: 1393–1512.

21. Fleiss JL, Levin B, Paik MC. *Statistical Methods for Rates and Proportions*. John Wiley & Sons Ltd, The Atrium, Southern Gate, Chichester, West Sussex, United Kingdom, https://onlinelibrary.wiley.com/doi/book/10.1002/0471445428 (2003).

22. Hemming K, Taljaard M. Key considerations for designing, conducting and analysing a cluster randomized trial. *Int J Epidemiol* 2023; 52: 1648–1658.

23. Li P, Redden DT. Comparing denominator degrees of freedom approximations for the generalized linear mixed model in analyzing binary outcome in small sample cluster-randomized trials. *BMC Medical Res Methodol* 2015; 15: 38.

24. Li P, Redden DT. Small sample performance of bias-corrected sandwich estimators for cluster-randomized trials with binary outcomes. *Stat Med* 2015; 34: 281–296.

25. Schluchter M, Elashoff J. Small-sample adjustments to tests with unbalanced repeated measures assuming several covariance structures. *J Stat Comput Sim* 1990; 37: 69–87.

26. Satterthwaite FE. An Approximate Distribution of Estimates of Variance Components. *Biometrics Bulletin* 1946; 2: 110–114.

27. Kenward MG, Roger JH. Small Sample Inference for Fixed Effects from Restricted Maximum Likelihood.
*Biometrics* 1997; 53: 983–997.

28. Fay MP, Graubard BI. Small-Sample Adjustments for Wald-Type Tests Using Sandwich Estimators.
*Biometrics* 2001; 57: 1198–1206.

29. Kauermann G, Carroll RJ. A Note on the Efficiency of Sandwich Covariance Matrix Estimation. *J Am Stat Assoc* 2001; 96: 1387–1396.

30. Morel J G., Bokossa M C., Neerchal N K. Small Sample Correction for the Variance of GEE Estimators.
*Biom* J 2003; 45: 395–409.

31. Mancl LA, DeRouen TA. A Covariance Estimator for GEE with Improved Small-Sample Properties.
*Biometrics* 2001; 57: 126–134.

32. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, https://cran.r-project.org/ (2023).

33. Sharples K, Breslow N. Regression analysis of correlated binary data: some small sample results for the estimating equation approach. *J Stat Comput Simul* 1992; 42: 1–20.

34. Le Vilain–Abraham F, Tavernier E, Dantan E, et al. Restricted mean survival time to estimate an intervention effect in a cluster randomized trial. *Stat Methods Med Res* 2023; 32: 2016–2032.

35. Landais M, Nay M-A, Auchable J, et al. Continued enteral nutrition until extubation compared with fasting before extubation in patients in the intensive care unit: an open-label, cluster-randomised, parallel-group, non-inferiority trial. *Lancet Respir Med* 2023; 11: 319–328.

36. Donner A, Klar N. *Design and Analysis of Cluster Randomization Trials in Health Research*. Wiley.
United Kingdom: John Wiley & Sons Ltd, The Atrium, Southern Gate, Chichester, West Sussex, United Kingdom, 2000.

37. Kahan BC, Blette BS, Harhay MO, et al. Demystifying estimands in cluster-randomised trials. *Stat Methods Med Res* 2024; 9622802241254197.

38. van der Laan M, Rubin D. Targeted Maximum Likelihood Learning. *Int J Biostat* 2006; 2(1).
39. Balzer LB, Zheng W, van der Laan MJ, et al. A new approach to hierarchical data analysis: Targeted maximum likelihood estimation for the causal effect of a cluster-level exposure. *Stat Methods Med Res* 2019; 28: 1761–1780.
# Appendices

## Appendix 1: Expected value of the individual-level covariate $X^{(m)}$

$$E(X_{ijl}^{(m)}) = E(A_{ijl}^{(m)} W_{ijl}^{(m)})$$

$$\square = E(A_{ijl}^{(m)})E(W_{ijl}^{(m)}),$$

$$\square \quad \square \quad \text{because } A_{ijl}^{(m)} \text{ and } W_{ijl}^{(m)} \text{ are independent}$$

$$\square \quad \square \quad \square$$

$$= \alpha_m E( (1 - R_{ijl}^{(m)}) G_{ijl}^{(m)} + R_{ijl}^{(m)} Y_{ijl} )$$

$$= \alpha_m [E(1 - R_{ijl}^{(m)})E(G_{ijl}^{(m)}) + E(R_{ijl}^{(m)})E(Y_{ijl})],$$

$$\square \quad \square \quad \text{because } Y_{ijl}, R_{ijl}^{(m)} \text{ and } G_{ijl}^{(m)} \text{ are independent}$$

$$\square \quad \square \quad \square$$

$$= \alpha_m [(1 - r_{i,m})p_i + r_{i,m}p_i]$$

$$= \frac{q_m}{p_i} p_i$$

$$= q_m$$

## Appendix 2: Variance of the individual-level covariate $X^{(m)}$

$$\text{Var}(X_{ijl}^{(m)}) = \text{Var}(A_{ijl}^{(m)} W_{ijl}^{(m)})$$

$$\square = \text{Var}(A_{ijl}^{(m)})\text{Var}(W_{ijl}^{(m)}) + \text{Var}(A_{ijl}^{(m)})E(W_{ijl}^{(m)})^2 + \text{Var}(W_{ijl}^{(m)})E(A_{ijl}^{(m)})^2$$

$$\square = \alpha_m(1 - \alpha_m)p_i(1 - p_i) + \alpha_m(1 - \alpha_m)p_i^2 + p_i(1 - p_i)\alpha_m^2$$

$$\square = \alpha_m p_i ((1 - \alpha_m)(1 - p_i) + (1 - \alpha_m)p_i + (1 - p_i)\alpha_m)$$

$$\square = \alpha_m p_i (1 - \alpha_m p_i)$$

$$\square = \frac{q_m}{p_i} p_i \left(1 - \frac{q_m}{p_i} p_i\right)$$

$$\square = q_m(1 - q_m)$$
Appendix 3: Correlation between the individual-level covariate $X^{(m)}$ and outcome $Y$

$$ \operatorname{Corr}(X_{ijl}^{(m)}, Y_{ijl}) = \frac{\operatorname{cov}(X_{ijl}^{(m)}, Y_{ijl})}{\sigma_{X_{ijl}^{(m)}} \sigma_{Y_{ijl}}} $$

$$ \square = \frac{\operatorname{cov}(A_{ijl}^{(m)} W_{ijl}^{(m)}, Y_{ijl})}{\sigma_{X_{ijl}^{(m)}} \sigma_{Y_{ijl}}} $$

$$ \square = \frac{\operatorname{cov}(A_{ijl}^{(m)} G_{ijl}^{(m)} - A_{ijl}^{(m)} R_{ijl}^{(m)} G_{ijl}^{(m)} + A_{ijl}^{(m)} R_{ijl}^{(m)} Y_{ijl}, Y_{ijl})}{\sigma_{X_{ijl}^{(m)}} \sigma_{Y_{ijl}}} $$

$$ \square = \frac{\operatorname{cov}(A_{ijl}^{(m)} R_{ijl}^{(m)} Y_{ijl}, Y_{ijl})}{\sigma_{X_{ijl}^{(m)}} \sigma_{Y_{ijl}}} $$

$$ \square = \frac{E(A_{ijl}^{(m)}) E(R_{ijl}^{(m)}) Var(Y_{ijl})}{\sigma_{X_{ijl}^{(m)}} \sigma_{Y_{ijl}}}, \quad \text{because } A_{ijl}^{(m)}, R_{ijl}^{(m)} \text{ and } Y_{ijl} \text{ are independent} $$

$$ \square = \frac{\alpha_m r_{i,m} p_i (1 - p_i)}{\sqrt{q_m (1 - q_m)} \sqrt{p_i (1 - p_i)}} $$

$$ \square = \frac{\alpha_m r_{i,m} \sqrt{p_i (1 - p_i)}}{\sqrt{q_m (1 - q_m)}} $$

$$ \square = \frac{\alpha_m \rho_m \frac{\sqrt{(1 - p_i)(1 - q_m)}}{\sqrt{\alpha_m (1 - p_i)}} \sqrt{p_i (1 - p_i)}}{\sqrt{q_m (1 - q_m)}} $$

$$ \square = \frac{\rho_m \sqrt{\alpha_m \sqrt{p_i}}}{\sqrt{q_m}} $$

$$ \square = \frac{\rho_m \sqrt{\frac{q_m}{p_i}} \sqrt{p_i}}{\sqrt{q_m}} $$

$$ \square = \rho_m $$
Appendix 4: Expression of the odds ratio associating the individual-level covariate $X_{ijl}^{(m)}$ and outcome $Y_{ijl}$

<table><thead><tr><th></th><th>Y<sub>ijl</sub> = 1</th><th>Y<sub>ijl</sub> = 0</th></tr></thead><tbody><tr><th>X<sub>ijl</sub><sup>(m)</sup> = 1</th><td>A</td><td>B</td></tr><tr><th>X<sub>ijl</sub><sup>(m)</sup> = 0</th><td>C</td><td>D</td></tr></tbody></table>

$$OR_m = \frac{A \ D}{C \ B} = \frac{\mathbb{P}(A)\mathbb{P}(D)}{\mathbb{P}(C)\mathbb{P}(B)}$$

We have:

To simplify the notations, we do not use curly brackets, so we write {$Y = 1$} as $Y = 1$ and simplified the index of the variable for example $Y_{ijl}$ as $Y$.

$$\mathbb{P}(A) = \mathbb{P}(Y=1 \cap X^{(m)}=1)$$

$$\square = \mathbb{P}(Y = 1 \cap A^{(m)} = 1 \cap R^{(m)} = 1) + \mathbb{P}(Y = 1 \cap A^{(m)} = 1 \cap G^{(m)} = 1 \cap R^{(m)} = 0)$$

$$\square = q_m r_{i,m} + p_i q_m - p_i q_m r_{i,m}$$

$$\mathbb{P}(B) = \mathbb{P}(Y=0 \cap X^{(m)}=1)$$

$$\square = \mathbb{P}(Y=0 \cap A^{(m)} = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 1)$$

$$\square = q_m - q_m r_m - p_i q_m + p_i q_m r_{i,m}$$

$$\mathbb{P}(C) = \mathbb{P}(Y=1 \cap X^{(m)}=0)$$

$$\square = \mathbb{P}(Y = 1 \cap A^{(m)} = 0) + \mathbb{P}(Y = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 0)$$

$$\square \square - \mathbb{P}(Y = 1 \cap A^{(m)} = 0 \cap R^{(m)} = 0 \cap G^{(m)} = 0)$$

$$\square = p_i - p_i q_m - q_m r_{i,m} + p_i q_m r_{i,m}$$

$$\mathbb{P}(D) = \mathbb{P}(Y=0 \cap X^{(m)}=0)$$

$$\square = \mathbb{P}(Y = 0 \cap A^{(m)} = 0) + \mathbb{P}(Y = 0 \cap R^{(m)} = 1)$$

$$\square = + \mathbb{P}(Y = 0 \cap G^{(m)} = 0) - \mathbb{P}(Y = 0 \cap R^{(m)} = 1 \cap G^{(m)} = 0)$$

$$\square \square - [\mathbb{P}(Y = 0 \cap A^{(m)} = 0 \cap R^{(m)} = 1) + \mathbb{P}(Y = 0 \cap A^{(m)} = 0 \cap G^{(m)} = 0)]$$

$$\square \square - \mathbb{P}(Y = 0 \cap A^{(m)} = 0 \cap R^{(m)} = 1 \cap G^{(m)} = 0)]$$

$$\square = 1 - q_m + q_m r_{i,m} - p_i + p_i q_m - p_i q_m r_{i,m}$$

With: $r_{i,m} = \rho_m \frac{\sqrt{1-q_m}}{\sqrt{\alpha_m(1-p_i)}}$

And so, we have:

$$OR_m = \frac{\mathbb{P}(A)\mathbb{P}(D)}{\mathbb{P}(C)\mathbb{P}(B)}$$

$$\square = \frac{(q_m r_{i,m} + p_i q_m - p_i q_m r_{i,m})(1 - q_m + q_m r_{i,m} - p_i + p_i q_m - p_i q_m r_{i,m})}{(q_m - q_m r_{i,m} - p_i q_m + p_i q_m r_{i,m})(p_i - p_i q_m - q_m r_{i,m} + p_i q_m r_{i,m})}$$
Appendix 5: Expression of the odds ratio associating the cluster-level covariate $X_{ij}^{(m)}$ and outcome $Y_{ijl}$

<table><thead><tr><th></th><th>Y<sub>ijl</sub> = 1</th><th>Y<sub>ijl</sub> = 0</th></tr></thead><tbody><tr><th>X<sub>ij</sub><sup>(m)</sup> = 1</th><td>A</td><td>B</td></tr><tr><th>X<sub>ij</sub><sup>(m)</sup> = 0</th><td>C</td><td>D</td></tr></tbody></table>

$$OR'_{m} = \frac{AD}{CB} = \frac{\mathbb{P}(A)\mathbb{P}(D)}{\mathbb{P}(C)\mathbb{P}(B)}$$

We have:

To simplify the notations, we do not use curly brackets, so we write {$Y = 1$} as $Y = 1$ and simplified the index of the variable for example $Y_{ijl}$ as $Y$.

$$\mathbb{P}(A) = \mathbb{P}(Y = 1 \cap X^{(m)} = 1)$$

$$\square = \mathbb{P}(U = 0 \cap V = 1 \cap A^{(m)} = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 1)$$

$$\square \square + \mathbb{P}(U = 0 \cap V = 1 \cap A^{(m)} = 1 \cap R^{(m)} = 1 \cap Z = 1)$$

$$\square \square + \mathbb{P}(U = 1 \cap Z = 1 \cap A^{(m)} = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 1)$$

$$\square \square + \mathbb{P}(U = 1 \cap Z = 1 \cap A^{(m)} = 1 \cap R^{(m)} = 1)$$

$$\square = p_i^2 \alpha + p_i \alpha r_m' \sqrt{\rho_{icc}} - p_i^2 \alpha r_m' \sqrt{\rho_{icc}}$$

$$\mathbb{P}(B) = \mathbb{P}(Y = 0 \cap X^{(m)} = 1)$$

$$\square = \mathbb{P}(U = 0 \cap V = 0 \cap A^{(m)} = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 1)$$

$$\square \square + \mathbb{P}(U = 0 \cap V = 0 \cap A^{(m)} = 1 \cap R^{(m)} = 1 \cap Z = 1)$$

$$\square \square + \mathbb{P}(= 1 \cap Z = 0 \cap A^{(m)} = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 1)$$

$$\square = p_i \alpha - p_i^2 \alpha - p_i \alpha r_m' \sqrt{\rho_{icc}} + p_i^2 \alpha r_m' \sqrt{\rho_{icc}}$$

$$\mathbb{P}(C) = \mathbb{P}(Y = 1 \cap X^{(m)} = 0)$$

$$\square = \mathbb{P}(U = 0 \cap V = 1 \cap A^{(m)} = 0) + \mathbb{P}(U = 1 \cap Z = 1 \cap A^{(m)} = 0)$$

$$\square \square + \mathbb{P}(U = 0 \cap V = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 0) + \mathbb{P}(U = 0 \cap V = 1 \cap R^{(m)} = 1 \cap Z = 0)$$

$$\square \square + \mathbb{P}(U = 1 \cap Z = 1 \cap R^{(m)} = 0 \cap G^{(m)} = 0)$$

$$\square \square - \mathbb{P}(U = 0 \cap V = 1 \cap A^{(m)} = 0 \cap R^{(m)} = 0 \cap G^{(m)} = 0)$$

$$\square \square - \mathbb{P}(U = 0 \cap V = 1 \cap A^{(m)} = 0 \cap R^{(m)} = 1 \cap Z = 0)$$

$$\square \square - \mathbb{P}(U = 1 \cap Z = 1 \cap A^{(m)} = 0 \cap R^{(m)} = 0 \cap G^{(m)} = 0)$$

$$\square = p_i - p_i^2 \alpha - p_i \alpha r_m' \sqrt{\rho_{icc}} + p_i^2 \alpha r_m' \sqrt{\rho_{icc}}$$

$$\mathbb{P}(D) = \mathbb{P}(Y = 0 \cap X^{(m)} = 0)$$

$$\square = \mathbb{P}(U = 0 \cap V = 0 \cap A^{(m)} = 0) + \mathbb{P}(U = 1 \cap Z = 0 \cap A^{(m)} = 0)$$

$$\square \square + \mathbb{P}(U = 0 \cap V = 0 \cap R^{(m)} = 0 \cap G^{(m)} = 0) + \mathbb{P}(U = 0 \cap V = 0 \cap R^{(m)} = 1 \cap Z = 0)$$

$$\square \square + \mathbb{P}(U = 1 \cap Z = 0 \cap R^{(m)} = 0 \cap G^{(m)} = 0) + \mathbb{P}(U = 1 \cap Z = 0 \cap R^{(m)} = 1)$$

$$\square \square + \mathbb{P}(U = 0 \cap V = 0 \cap A^{(m)} = 0 \cap R^{(m)} = 0 \cap G^{(m)} = 0)$$

$$\square \square - \mathbb{P}(U = 0 \cap V = 0 \cap A^{(m)} = 0 \cap R^{(m)} = 1 \cap Z = 0)$$

$$\square \square - \mathbb{P}(U = 1 \cap Z = 0 \cap A^{(m)} = 0 \cap R^{(m)} = 0 \cap G^{(m)} = 0)$$

$$\square \square - \mathbb{P}(U = 1 \cap Z = 0 \cap A^{(m)} = 0 \cap R^{(m)} = 1)$$

$$\square = 1 - p_i \alpha + p_i \alpha r_m' \sqrt{\rho_{icc}} - p_i + p_i^2 \alpha - p_i^2 \alpha r_m' \sqrt{\rho_{icc}}$$
And so, we have,

$$OR'_{m} = \frac{\mathbb{P}(A)\mathbb{P}(D)}{\mathbb{P}(C)\mathbb{P}(B)}$$

$$\Box = \frac{(p_i q'_m + q'_m r'_m \sqrt{\rho} - p q'_m r'_m \sqrt{\rho})(1 - q'_m + q'_m r'_m \sqrt{\rho} - p_i + p_i q'_m - p_i q'_m r'_m \sqrt{\rho})}{(q'_m - q'_m r'_m \sqrt{\rho} - p_i q'_m + p_i q'_m r'_m \sqrt{\rho})(p_i - p_i q'_m - q'_m r'_m \sqrt{\rho} + p_i q'_m r'_m \sqrt{\rho})}$$
# List of Tables and Figures

## Table 1: Simulation plan parameters

## Table 2: Performance measures

**Figure 1:** Proportion of scenarios with type I error rate considered liberal (type I error rate > 7%, dark blue), conservative (type I error rate < 3%, light blue) or acceptable (type I error rate 3% to 7%, blue) for scenarios with no treatment effect. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.

**Figure 2:** Bias in scenarios with no treatment effect. Results are displayed according to the number of clusters per arm (columns), intra-cluster correlation coefficient (rows) and total number of covariates (dot color).

**Figure 3:** Proportion of scenarios with relative error of the standard error (SE) rate considered over-estimated (relative error of the SE rate > 10%, dark blue), under-estimated (relative error of the SE rate < -10%, light blue) or acceptable (relative error of the SE rate -10% to 10%, blue) for scenarios with no treatment effect. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.

**Figure 4:** Proportion of scenarios with coverage rate considered liberal (coverage rate < 93%, dark blue), conservative (coverage rate > 97%, light blue) or acceptable (coverage rate 93% to 97%, blue) for scenarios with non-null treatment effects. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.

**Figure 5:** Bias in scenarios with non-null treatment effect. Results are displayed according to the number of clusters per arm (columns), intra-cluster correlation coefficient (rows) and total number of covariates (dot color).

**Figure 6:** Proportion of scenarios with relative error of the SE rate considered over-estimated (relative error of the SE rate > 10%, dark blue), under-estimated (relative error of the SE rate < -10%, light blue) or acceptable (relative error of the SE rate -10% to 10%, blue) for scenarios with non-null treatment effect. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.

**Figure 7:** Analysis of the Ambroisie non-inferiority cluster randomized trial (non-inferiority margin of 0.1). The original result (Ambroisie analysis), derived using a variance correction is compared to GEE and GLMM approaches using Gaussian, binomial and Poisson distribution without (Crude) and with (Adjusted) adjustment on binary covariates: obesity status, use of corticoids, chronic respiratory status, heart disease status and university status of the hospital.
<table><caption>Table 1: Simulation plan parameters</caption><thead><tr><th>Parameter</th><th>Number</th><th>Notation</th><th>Value</th></tr></thead><tbody><tr><td>Pair of prevalences</td><td>6</td><td>(p<sub>1</sub>, p<sub>0</sub>)</td><td>No treatment effect:<br>(0.5, 0.5)<br>(0.2, 0.2)<br>(0.1, 0.1)<br>Non-null treatment effect:<br>(0.5, 0.4)<br>(0.2, 0.15)<br>(0.1, 0.05)</td></tr><tr><td>Number of clusters per arm</td><td>5</td><td>k</td><td>5<br>10<br>20<br>30<br>40</td></tr><tr><td>Mean cluster size</td><td>9</td><td>n̄</td><td>25, 280 when k = 5<br>15, 200 when k = 10<br>10, 150 when k = 20, 30<br>50 when k = 40</td></tr><tr><td>Intra-cluster correlation coefficient</td><td>3</td><td>ρ</td><td>0.001<br>0.01<br>0.05</td></tr><tr><td>Number of individual-level covariates,</td><td>3</td><td>S</td><td>0, 1, 3<br>Scenarios with 1 covariate:<br>OR<sub>1</sub> = 2.5, β<sub>(1),1</sub> = 2.5, β<sub>(1),0</sub> = 1<br>Scenarios with 3 covariates:<br>OR<sub>1</sub> = 2.5, OR<sub>2</sub> = 1.5, OR<sub>3</sub> = 1,<br>β<sub>(1),1</sub> = 2.5, β<sub>(2),1</sub> = 2.5, β<sub>(3),1</sub> = 2.5<br>β<sub>(1),0</sub> = 1, β<sub>(2),0</sub> = 1, β<sub>(3),0</sub> = 1</td></tr><tr><td>Number of cluster-level covariates</td><td>3</td><td>S'</td><td>0, 1, 2<br>Scenarios with 1 covariate:<br>{OR<sub>1</sub> = 1.2, if ρ = 0.05/0.01<br>OR<sub>1</sub> = 1.1, if ρ = 0.001<br>Scenarios with 2 covariates:<br>{OR<sub>1</sub> = 1.2, OR<sub>2</sub> = 1, if ρ = 0.05/0.01<br>OR<sub>1</sub> = 1.1, OR<sub>2</sub> = 1, if ρ = 0.001</td></tr></tbody></table>

OR, odds ratio

$β_{(m),i}$ the parameter quantifying the association between the $m^{th}$ individual-level covariate with the inclusion status in the $i^{th}$ arm
<table><caption>Table 2: Performance measures</caption><thead><tr><th>Performance measures</th><th>Definition</th><th>Estimation</th></tr></thead><tbody><tr><td>Type I error*</td><td>Pr(θ &notin; [θ<sub>low</sub>; θ<sub>upp</sub>] | θ = 0)</td><td>1/n<sub>sim</sub> &sum;<sub>b=1</sub><sup>n<sub>sim</sub></sup> 1 (θ &notin; [θ<sub>low</sub>; θ<sub>upp</sub>])</td></tr><tr><td>Coverage rate*</td><td>Pr(θ̂<sub>low</sub> &le; θ &le; θ̂<sub>upp</sub>)</td><td>1/n<sub>sim</sub> &sum;<sub>b=1</sub><sup>n<sub>sim</sub></sup> 1 (θ̂<sub>low,b</sub> &le; θ &le; θ̂<sub>upp,b</sub>)</td></tr><tr><td>Bias</td><td>E[θ̂] - θ</td><td>1/n<sub>sim</sub> &sum;<sub>b=1</sub><sup>n<sub>sim</sub></sup> θ̂<sub>b</sub> - θ</td></tr><tr><td>Standard error relative bias (%)**</td><td>100 (ModSE / EmpSE - 1)</td><td>100 (ModSE / EmpSE - 1)</td></tr></tbody></table>

$$ \text{\*\*ModSE: } \sqrt{E[\widehat{Var}(\hat{\theta})]} \text{ and } \widehat{ModSE} = \sqrt{\frac{1}{n_{sim}} \sum_{b=1}^{n_{sim}} \widehat{Var}(\hat{\theta}_b)} $$

$$ \text{\*\*EmpSE: } \sqrt{\widehat{Var}(\hat{\theta})} \text{ and } \widehat{EmpSE} = \sqrt{\frac{1}{n_{sim}-1} \sum_{b=1}^{n_{sim}} (\hat{\theta}_b - \bar{\theta})^2} $$

<footnote>* The confidence intervals of the b$^{th}$ iteration were calculated as: $CI_b = [\hat{\theta}_{low,b} ; \hat{\theta}_{upp,b}]$, with $\hat{\theta}_{low,b} = \hat{\theta}_b - t_{\frac{\alpha}{2},DF}(\hat{SE}_b)$ and $\hat{\theta}_{upp,b} = \hat{\theta}_b + t_{\frac{\alpha}{2},DF}(\hat{SE}_b)$ and DF is the degree of freedom as determined in section “Standard error”.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Proportion of scenarios with type I error rate considered liberal (type I error rate > 7%, dark blue), conservative (type I error rate < 3%, light blue) or acceptable (type I error rate 3% to 7%, blue) for scenarios with no treatment effect. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Bias in scenarios with no treatment effect. Results are displayed according to the number of clusters per arm (columns), intra-cluster correlation coefficient (rows) and total number of covariates (dot color).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3: Proportion of scenarios with relative error of the standard error (SE) rate considered over-estimated (relative error of the SE rate > 10%, dark blue), under-estimated (relative error of the SE rate < -10%, light blue) or acceptable (relative error of the SE rate -10% to 10%, blue) for scenarios with no treatment effect. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Proportion of scenarios with coverage rate considered liberal (coverage rate < 93%, dark blue), conservative (coverage rate > 97%, light blue) or acceptable (coverage rate 93% to 97%, blue) for scenarios with non-null treatment effects. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Bias in scenarios with non-null treatment effect. Results are displayed according to the number of clusters per arm (columns), intra-cluster correlation coefficient (rows) and total number of covariates (dot color).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Proportion of scenarios with relative error of the SE rate considered over-estimated (relative error of the SE rate > 10%, dark blue), under-estimated (relative error of the SE rate < -10%, light blue) or acceptable (relative error of the SE rate -10% to 10%, blue) for scenarios with non-null treatment effect. Each row is associated with one simulation parameter: number of clusters per arm, intra-cluster correlation coefficient, total number of covariates, pairs of prevalence, mean cluster size.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Analysis of the Ambroisie non-inferiority cluster randomized trial (non-inferiority margin of 0.1). The original result (Ambroisie analysis), derived using a variance correction is compared to GEE and GLMM approaches using Gaussian, binomial and Poisson distribution without (Crude) and with (Adjusted) adjustment on binary covariates: obesity status, use of corticoids, chronic respiratory status, heart disease status and university status of the hospital.</figcaption></figure>