The ICCs have been estimated from the data that was generated in the Targeted Malaria Elimination study (TME) with mass drug administrations (MDA) on falciparum malaria in South-East Asia [12]. Following vector control activities, community-based case management and intensive community engagement, restricted randomization was conducted within village pairs to select 8 villages to receive early MDA and 8 villages as controls. After 12 months the control villages received deferred MDA. The MDA comprised 3 monthly rounds of 3 daily doses of dihydroartemisinin–piperaquine and, except in Cambodia, a single low dose of primaquine. Cross-sectional surveys of the entire population of each village at quarterly intervals using ultrasensitive quantitative PCR (uPCR) were used to detect *Plasmodium* infections. The overall aim of the study was to assess the duration of effectiveness of MDA on falciparum parasitaemia incidence and prevalence in 16 remote village populations, 4 villages each in Myanmar, Vietnam, Cambodia and Laos. The sample size, 4 village clusters per country, was chosen mainly for operational and practical reasons. The detailed methods of the TME study have been published [12,13,14,15,16].

### Definitions of outcomes

Defining a new *P. vivax* infection from longitudinally collected data is more complicated than *P. falciparum* as *P. vivax* infections recur frequently. Recurrences of *P. vivax* infections after treatment of the blood stage infection can be due to recrudescence, relapse or reinfection. The ICCs are based on baseline prevalence of *P. vivax*/*falciparum* infections and the cumulative incidence of detected *P. vivax*/*falciparum* parasitaemias at each quarterly surveys based on uPCR results.

### Estimation of *P. vivax*/*falciparum* prevalence and incidence over a 12-month period

The cumulative incidence of *P. vivax/falciparum* parasitaemias over the 12-month period was calculated based on uPCR results collected at month 0, 3, 6, 9 and 12. A participant was considered to have a recurrent *P. vivax* infection if there were two or more positive uPCR results during the 12 months follow up period. As consecutive positive uPCR tests could be due to a re-infection following a new mosquito bite or a continuous infection which is likely to be due to persistence in *P. falciparum* infections and a relapse in *P. vivax* infections. To address this uncertainty, an “episode” was defined in two ways. In the first approach each positive uPCR test was considered as a separate episode (i.e. reinfections). In the second, consecutive positive uPCR results were considered to belong to the same continuous infection (persistent or relapsing asymptomatic parasitaemia). The ICCs for the second approach are presented in the Additional file 1.

### Statistical methodology

The outcomes of interest for the ICC estimation are the prevalence and the incidence of *P. falciparum* and *P. vivax* infections. A logistic regression model is the most relevant model for prevalence while the Poison model is the most natural model when modeling incidence as the outcome of interest. The basic ICC formula presented above (\(\rho = \frac{{\sigma_{b}^{2} }}{{\sigma_{b}^{2} + \sigma_{w}^{2} }}\)) refers to a case where two-level hierarchical data are of interest. In practice, often several levels of clustering are available and of interest. The hierarchical structure of the data in the TME study included 4 levels: longitudinal data on infection status (level 1) collected repeatedly for each individual (level 2) who belonged to a village (level 3) which was located in a country (level 4). However, country specific ICCs were estimated because there was considerable heterogeneity in baseline *P. falciparum*/*P. vivax* prevalence between countries. In this case, the level of country is not considered. The model for estimating ICC for prevalence is reduced to 2 levels in each country because ICCs were estimate at baseline only and, therefore, each individual contributes only one observation at baseline. By contrast, for the estimation of ICC for incidence, multiple outcomes were aggregated, i.e. each individual had one observation for the outcome counts over time and the exposure time was aggregated for each individual. Two-level hierarchical models were fitted to estimate country specific ICCs for both prevalence and incidence with a village as unit of randomization.

Methodological approaches have been developed describing procedures used to compute ICC values applicable to models with multiple hierarchical levels that include logistic regression models as well as other generalized linear models such as the Poisson regression models [8, 10, 11]. The ICC values can be estimated from model equations as exact estimation methods or through use of simulations. A latent variable approach is another method for estimating ICC from logistic (logit link-scale) regression models. The link-scale is often considered to be of interest for prevalence, because the estimates of the individual outcomes are performed on the underlying latent scale [17, 18]. Nakagawa et al. provided comprehensive methods for calculating ICCs using both exact and latent variable methods for logistic and Poisson models [18]. However, Austin et al. have shown that there is no latent response formulation for Poisson models and such a model is, therefore, not included in this paper. The exact estimation method and the simulation method were utilized in the estimation of the ICC values for both incidence (Poisson model) and prevalence (logit model) of *P. falciparum* and *P. vivax*. In addition, the latent variable estimates of ICC for logit model are provided in the additional material for the estimation of ICCs for prevalence of *P. falciparum* and *P. vivax*. ICC values for the model are provided with and without the covariates sex and age because they were independently associated with the outcome [12]. The estimation of the country specific ICC is the main focus of this article. However, the overall ICC is also included for prevalence of *P*. *falciparum* and *P. vivax* using the latent variable method. For prevalence ICC uses the baseline prevalence as this is the time-point measure most often used in sample size calculations. The statistical methodology for estimating ICC from a random effects logistic model using the exact calculation, simulation-based and latent variable method is introduced. In addition the methodology for estimating ICCs from the random effect Poisson model using exact calculation and simulation-based methods is described [8].

For prevalence outcome, consider a logistic model for the outcome \(Y_{ij}\), where \(i\) denotes an individual and \(j\) denotes a cluster then:

$$Y_{ij} \sim {\text{Bernoulli}}\left( {p_{ij} } \right)$$

(1)

where \(p_{ij}\) is the probability of experiencing an outcome for individual \(i\) in cluster \(j\). And the logistic regression model is fitted as a generalized linear mixed model (GLMM) with logit link function as:

$${\text{logit (}}p_{ij} ) {\text{ = log}}\left( {\frac{{p_{ij} }}{{1 - p_{ij} }}} \right) = \alpha_{j} + \beta X_{ij}$$

(2)

where \(X_{ij}\) refers to covariates such as age and sex measured on the individual \(i\) in cluster \(j\) and \(\alpha_{j}\) is a cluster-specific random effect such that \(\alpha_{j}\) follows a normal distribution with a mean of 0 and variance equal to \(\sigma_{\alpha }^{2}\).

The exact ICC for prevalence using logistic model is calculated as follows [10]

$${\text{ICC}} = \frac{{\left[ {\sigma_{\alpha }^{2} p_{ij}^{2} /\left( {1 + \exp \left( {2\beta X + \sigma_{\alpha }^{2} } \right)} \right)} \right]}}{{\left[ {\sigma_{\alpha }^{2} p_{ij}^{2} /\left( {1 + \exp \left( {2\beta X + \sigma_{\alpha }^{2} } \right)} \right)} \right] + \left[ {p_{ij} \left( {1 - p_{ij} } \right)} \right]}} ,{\text{ where}}\;p_{ij} = \frac{{\exp \left( {\beta X} \right)}}{{1 + \exp \left( {\beta X} \right)}}$$

(3)

where \(\sigma_{\alpha }^{2}\) (the random effect variance for the level of interest) and \(\beta\) [the log (odds ratio)] are estimated from the model and \(X\) refers to covariates such as age and sex.

The calculation of ICC as a postestimation estimate from software is provided using the latent variable approach for logistic model. In the logistic model, the underlying logistic error distribution has a constant variance \(\frac{{\pi^{2} }}{3}\) which was used as residual variance when calculating ICC from a logistic model and then the latent ICC is given by \(ICC = \frac{{\sigma_{\alpha }^{2} }}{{\sigma_{\alpha }^{2} + \frac{{\pi^{2} }}{3}}}\), where \(\sigma_{\alpha }^{2}\) is the between cluster variance for the binary outcome. Goldstein et al. [11] suggested that the latent variable approach to estimate the ICC is only appropriate when the binary outcome can be an underlying continuous latent variable. However, this is the version of ICC that is readily obtained in most software including Stata. In a Poisson model, the error variance is not constant and depends on the covariates included in the model. Consider a Poisson model for the outcome \(Y_{ij}\), where \(i\) denotes an individual and \(j\) denotes a cluster then:

$$Y_{ij} \sim {\text{Poisson}}\left( {\lambda_{ij} } \right)$$

(4)

And the Poisson regression model is fitted as a generalized linear mixed model (GLMM) with log link function as:

$${ \log }\left( {\lambda_{ij} } \right) = \alpha_{j} + \beta X_{ij}$$

(5)

where \(X_{ij}\) refers to covariates such as age and sex measured on the individual \(i\) in cluster \(j\); \(\lambda_{ij}\) is an estimate of the expected number of outcome events for individual \(i\) in cluster \(j\); and \(\alpha_{j}\) is a cluster-specific random effect such that \(\alpha_{j}\) follows a normal distribution with a mean of 0 and variance equal to \(\sigma_{\alpha }^{2}\). The estimate of the ICC from exact calculation is calculated as follows [8]:

$${\text{ICC}} = \frac{{\left[ {\exp \left( {2\beta X + 2\sigma_{\alpha }^{2} } \right) - \exp \left( {2\beta X + \sigma_{\alpha }^{2} } \right)} \right]}}{{\left[ {\exp \left( {2\beta X + 2\sigma_{\alpha }^{2} } \right) - \exp \left( {2\beta X + \sigma_{\alpha }^{2} } \right)} \right] + \left[ {\exp \left( {\beta X + \frac{{\sigma_{\alpha }^{2} }}{2}} \right)} \right]}}$$

(6)

where \(\sigma_{\alpha }^{2}\) (the random effect variance for the level of interest) and \(\beta\)(the log (incidence rate ratio)) are estimated from the model and \(X\) refers to covariates such as age and sex. The simulation procedures are detailed in Austin et al. [8].

The simulation-based algorithm for both prevalence and incidence proceeds as follows:

- 1.
Fit a multilevel logistic or Poisson model to an existing dataset. If the dataset is not available, one can use desired parameters estimated from previous studies to generate a dataset that mimics the original study data.

- 2.
Simulate a large number say, *M*, cluster-level random effects (\(\alpha_{m} , m = 1, \ldots ,M)\) from a mixed effect model distribution obtained from the model fitted in step 1. Typically *M *= 1000; 5000; or 10, 000 simulations but higher numbers of *M* are also possible. (5000 simulations or more may take days, weeks or months depending on the data and the fitted model.)

- 3.
Use each of the simulated mixed effects model drawn in step 2 to compute the predicted means denoted as: (\(\left( {\hat{E}\left( {Y_{m} } \right)} \right)\), and variances denoted as: \(\hat{V}\left( {Y_{m} } \right)\),

- 4.
The predicted mean is calculated as: \(\hat{E}\left( {Y_{m} } \right) = \hat{p}_{m} = \frac{{\exp \left( {\hat{\beta }X + \alpha_{m} } \right)}}{{1 + \exp \left( {\hat{\beta }X + \alpha_{m} } \right)}}\) and the variance is \(\hat{V}\left( {Y_{m} } \right) = \hat{p}_{m} \left( {1 - \hat{p}_{m} } \right)\) for logistic, where \(\hat{p}_{m}\) is the probability of outcome predicted from the random effects model.

- 5.
The predicted mean and variance for the Poisson model are calculated as: \(\hat{E}\left( {Y_{m} } \right) = \hat{V}\left( {Y_{m} } \right) = \hat{\lambda }_{m} = \exp \left( {\hat{\beta }X + \alpha_{m} } \right)\).

- 6.
Then compute the between cluster variance, \(\hat{\sigma }_{\alpha }^{2} = V\left( {\hat{E}\left( {Y_{m} } \right)} \right)\), variance in the distribution of \(\left( {\hat{E}\left( {Y_{m} } \right)} \right)\) from the mean of the \(\left( {\hat{E}\left( {Y_{m} } \right)} \right)\) from the M simulations.

- 7.
Calculate the error variance as the mean of the estimated outcome variances \(E\left( {\hat{V}\left( {Y_{m} } \right)} \right) = \frac{1}{M}\sum\nolimits_{m = 1}^{M} {\hat{V}\left( {Y_{m} } \right)}\)

- 8.
Thus, the between cluster variance is \(V\left( {\hat{E}\left( {Y_{m} } \right)} \right)\) while the total variance is \(V\left( {\hat{E}\left( {Y_{m} } \right)} \right) + \frac{1}{\text{M}}\sum\nolimits_{m = 1}^{M} {\hat{V}\left( {Y_{m} } \right)}\).

- 9.
The estimate of the ICC from simulation-based approach is calculated as follows [8]:

$${\text{ICC}} = \frac{{V\left( {\hat{E}\left( {Y_{m} } \right)} \right)}}{{V\left( {\hat{E}\left( {Y_{m} } \right)} \right) + \frac{1}{\text{M}}\sum\limits_{m = 1}^{M} {\hat{V}\left( {Y_{m} } \right)} }}$$

The confidence intervals for the ICCs were calculated by using bootstrapped samples to estimate the standard error. All analyses including simulations and bootstrapping were performed in Stata 15.