{'fraction_non_alphanumeric': 0.035484220213284255, 'fraction_numerical': 0.011963406052076003, 'mean_word_length': 4.550931490384615, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Frequentist inference for cluster randomised trials with multiple primary outcomes

The use of a single primary outcome is generally either recommended or required by many influential randomised trial guidelines to avoid the problem of “multiple testing”. Without correction, the probability of rejecting at least one of a set of null hypotheses (the family-wise error rate) is often much greater than the nominal rate of any single test so that statistics like p-values and confidence intervals have no reliable interpretation. Cluster randomised trials though may require multiple outcomes to adequately describe the effects of often complex and multi-faceted interventions. We propose a method for inference for cluster randomised trials with multiple outcomes that ensures a nominal family-wise error rate and produces simultaneous confidence intervals with nominal “family-wise” coverage. We adapt the resampling-based stepdown procedure of Romano and Wolf (2005) using a randomisation-test approach within a generalised linear model framework. We then adapt the Robbins-Monro search procedure for confidence interval limits proposed by Garthwaite and Buckland (1996) to this stepdown process to produce a set of confidence intervals. We show that this procedure has nominal error rates and coverage in a simulation-based study of parallel and stepped-wedge cluster randomised studies and compare results from the analysis of a real-world stepped-wedge trial under both the proposed and more standard analyses.


Introduction
The requirement to state a single primary outcome for a randomised controlled trial has become accepted, even required, practice. For example, the influential CONSORT statement on clinical trials requires the pre-specification of a primary outcome, which they describe as "the pre-specified outcome considered to be of greatest importance to relevant stakeholders", and recommends against multiple primary outcomes [The Lancet, 2010]. The reason for this is to ensure appropriate control of the 'false discovery rate' when using null hypothesis significance testing [Wason et al., 2014]. If there are multiple outcomes each with their own associated treatment effect being tested separately, then we are implicitly testing a family of null hypotheses against an alternative that at least one of them is false.
Without correction, the type I error rate for this family of null hypotheses will be much greater than the nominal rate of any single test [Romano and Wolf, 2005a]. Indeed, the CONSORT statement notes that that multiple primary outcomes are not recommended as it "incurs the problem of multiplicity of analyses" [The Lancet, 2010].
Cluster randomised trials are a widely used method to evaluate interventions applied to groups of people, such as clinics, schools, or villages. Often these interventions target 'higher level' processes and can be complex in nature [Murray, 1998, Hemming et al., 2015, Eldridge and Kerry, 2012. Recent examples from our own work include an incentive scheme to improve implementation of a broad package of education and activities designed to improve employee health in the workplace [Thrive at Work Wellbeing Programme Collaboration, 2019a], or a community health worker programme targeting multiple health conditions [Dunbar et al., 2018]. The effects of such complex interventions cannot be adequately summarised by a single outcome. 'Fudging' this process by creating a composite outcome is undesirable since it requires applications of arbitrary weights across outcomes and discards information by collapsing a multivariate outcome to a univariate one. The requirement for a single primary outcome therefore clashes with the needs of many cluster randomised trials. The solution is to ensure appropriate methods are used where there are multiple outcomes of interest rather than restricting the outcomes from which we can make inferences. However, the question of appropriate analysis for randomised trials, and particularly cluster randomised trials, with multiple outcomes can be contentious and complex.
The Food and Drug Administration (FDA), the main regulatory body for medicines in the United States, declares that "If the purpose of the trial is to demonstrate effects on all of the designated primary variables, then there is no need for adjustment of the type I error" [US Food and Drug Administration, 2017]. They also identify a "gatekeeping" approach where "statistical significance" on a primary outcome is required before a second one can be analysed and state this does not need correction for multiple testing. However, in both cases, the type I error for the two (or more) outcomes jointly is unlikely to be nominal with this procedure and their validity is debatable. Other authors differentiate aiming to declare "statistical significance" on at least one of a group of null hypotheses to declaring it for all, and propose different solutions for both [Lafaye de Micheaux et al., 2014, Rubin, 2021.
Where a correction for multiple testing is deemed neccesary, we can divide solutions into: (i) multivariate methods that model the joint distribution of the outcomes, which is particularly favoured by Bayesian practitioners [Gelman et al., 2012]; and (ii) univariate solutions that aim to ensure inferential statistics for a set of estimands collectively have the appropriate Frequentist properties [Farcomeni, 2008]. Despite the different approaches and guidance, Wason et al. [2014] estimated that only around half of all randomised trials with multiple outcomes or arms corrected for multiple testing and outcomes. No evidence is available on the use of corrections for multiple testing in cluster randomised trials specif-ically, but there are few, if any, comprehensive discussions of methods in this area currently available. No guidance is available on methods for correcting confidence intervals or if it is needed; the FDA note that it is complex and beyond the scope of their advice.
This article adapts and evaluates a univariate method for the use of multiple primary outcomes in cluster randomised trials. We do not aim to resolve the debate about the necessity of these methods, rather to provide a novel method that does ensure a nominal type I error rate and confidence set coverage when multiple outcomes are used in cluster trials, should one be required. We do not consider the problem with the goal of declaring "statistical significance", although we do use terms such as "error rate" for convenience.
Our contention is that even in the absence of declarations of "statistical significance", statistics such as p-values and confidence intervals must possess the appropriate Frequentist properties in order to be interpretable. 95% confidence intervals are presented in applied research on the basis that they have a coverage of 95%, if instead they had unknown coverage that might range, say, between 50% and 99% then they would have little ability to support principled inference from statistical analysis of randomised trials and other studies. If a p-value of 0.03 could mean the null model(s) would produce datasets with more extreme test statistics in anywhere from 2% to 20% of repetitions, then it is similarly unusable. Our aim therefore is to derive p-values and confidence intervals with a nominal type I error rate for a family of null hypotheses and coverage of confidence sets for the parameters of interest.
2 Family wise error rates

The multiple testing problem
We first suppose that data X are generated from some probability distribution P , which belongs to some family of probability distributions Ω. The family Ω could be a parametric, semi-parametric, or non-parametric model. The multiple testing problem arises when we have a set of hypotheses H j versus H j for j = 1, ..., J, following the notation of Romano and Wolf [2005a]. Each of these hypotheses is a subset ω j ⊂ Ω and is equivalent to testing P ∈ ω j against P ∈ ω j . So for any subset K ⊂ 1, ..., J, H K = ∩ j∈K H j is the hypothesis that P ∈ ∩ j∈K ω j . We assume each null hypothesis H j is based on a test statistic T j ; we denote the α-quantile of the distribution of T j as c j (α, P ). In a traditional null hypothesis testing framework we "reject" H j in favour of H j at the α level, if T j ≥ c j (1 − α, P ), which clearly has probability α. Conversely, the p-value p j of the test is where T j = c j (1 − p j , P ), so that the probability of observing P r(T j > c j (1 − p j , P )|H j ) = p j . The family-wise error rate (FWER) of this set of hypotheses is the probability of "rejecting" at least one true null hypothesis. That is, if I = I(P ) ⊂ 1, ..., J are the indices of the true null hypotheses, so j ∈ I if and only if P ∈ ω j , then the FWER is the probability under P of rejecting any H j∈I , i.e. P r(∪ j∈I T j > c j (1 − α, P )), which should be α.

Solutions
Solutions to the multiple testing problem aim to ensure that F W ER ≤ α. Control over the FWER is said to be strong if it holds for any combination of true and false null hypotheses, and weak if it only holds when all null hypotheses are true [Dudoit et al., 2003]. Several approaches exist to control the FWER. The Bonferroni method is probably most well known, which sets the critical value for the test of the null hypothesis to be c j (1 − α/J, P ). However, while this method exerts strong control over the FWER. It is very conservative and assumes no correlation between the test statistics. Holm [1979] proposed a less conservative 'stepdown' approach to multiple testing. One orders the test statistics from largest to smallest and then compares the largest statistic to the critical value c j (1 − α/J, P ). If the test statistic is larger than this value, then the null hypothesis is rejected, otherwise we do not reject any null hypothesis and stop. If we rejected, then the next largest test statistic is compared to c j (1−α/(J −1), P ), and again it is either rejected, or we do not reject all remaining null hypotheses and stop, and so forth. The stepdown method is less conservative than the Bonferroni method, but it can still be conservative as it does not take into account the dependence structure in the data and of the test statistics. Romano and Wolf [2005a,b] developed a resampling based version of the stepdown method. Being based on randomisation or permutation tests provides an advantage as it allows for arbitrary dependence structures in the data while maintaining strong control of the FWER. Their process is optimal in a maximin sense. These properties, and the use of randomisation tests, are useful for cluster randomised trials with multiple outcomes as the construction of appropriate test statistics can be difficult owing to the difficulty of estimating reliable standard errors with the potentially complex covariance structure and often small number of clusters [Watson et al., 2021, Leyrat et al., 2018. Furthermore, as we discuss in Section 4, this procedure provides a means to construct confidence intervals.
Other methods for controlling the FWER have been proposed, which are surveyed elsewhere [Farcomeni, 2008].
We describe the general stepdown procedure of Romano and Wolf firstly in terms of accepting or rejecting each null hypothesis at an α-level. We let c K (α, P ) denote an α-quantile of the distribution of the statistic: for any subset of null hypotheses K. We also denote T |r| as the rth largest test statistic so that corresponding to hypotheses H |1| , H |2| , ..., H |J| . Then the idealised algorithm is: . . .
In this procedure it is assumed the critical values are known. However, as Romano and Wolf show, they can be replaced by estimates from randomisation tests. In particular, the null hypothesis implies that the distribution of X is invariant under a set of transformations in A, which has L elements. So, aX and X have the same distribution for all a ∈ A whenever X has distribution P ∈ ω.
Our observed test statistics with our sample data are T j (X). The test statistic generated by the lth permutation is T j (a l X) for a l ∈ A and l = 1, ..., L. For each permutation we can determine the test statistic as in Equation (1) as T K,l = max j∈K T j (a l X). As before we denote T K,|r| as the rth largest of all the permutational test statistics {T K,l ; l = 1, ..., L} of the test statistics from the permutations. Then our estimator for the critical value is: We can see how this procedure produces p-values for a two-sided hypothesis that also maintains the FWER for a given α [Romano and Wolf, 2016], in particular: For a one-sided test we would not use the absolute values of the test statistics.
Often the size of A can be very large, and increases exponentially with the number of clusters. A Monte Carlo approach can be used that instead generates a random subset of A of fixed sized in order to generate realisations of the test statistics. If we conduct M such permutations then the estimator of the p-value for a given null hypothesis versus some alternative isp Obtaining p-values in this way is described in detail by Romano and Wolf [2016]. In the next section we discuss permuatation tests in the context of cluster randomised trials.

Permutation tests and p-values for cluster trials
We next introduce a generalised linear mixed model commonly used in the analysis of cluster randomised trials (e.g. Hooper et al. [2016]). We denote Y ict as the outcome of the ith individual, i = 1, ..., N , in cluster c = 1, ..., C at time t = 1, ..., T . We do not restrict the outcome, it could be continuous or discrete. We also specify the linear predictor: where D ct is an indicator for whether cluster c has received the intervention at time t and so δ is the parameter of interest, our "treatment effect". We also have a vector of individual and/or cluster-level covariates, X ict , which may also contain temporal fixed effects. The parameter θ ct represents a general 'random-effect' term that captures the within cluster and cluster-time correlation, although we do not provide a specific structure here. The overall model is then where h(.) is a link function. For example, P could be a Binomial distribution and h(.) the logistic link function.
Gail et al. [1996] provided the first extensive examination of permutation tests for cluster-based study designs. Their work principally used unweighted differences of cluster means as the basis of permutation tests. In particular, ifȲ .ct = i Y ict then one can base inference on the mean difference: Under the null hypothesis that H 0 : U = 0 the clusters are exchangeable, so we can re-assign them to a new allocation a ∈ A. In this way critical values of the distribution of U under this null, or associated p-values, can be generated. This statistic is the basis of a statistical package for inference from stepped-wedge cluster randomised trials using randomisation tests (cptest; Gallis et al. [2018]). The authors of that package also implement a "withinperiod" version of the analysis which takes estimates "within period" treatment effects as: and then takes a weighted average of them based on the number of clusters in treatment and control conditions in each time period. One can extend these ideas to allow for covariate adjustment [Gail et al., 1996]. Braun and Feng [2001] examine optimal permutation tests for cluster randomised trials.
The test statistics they consider are weighted sums of residuals. In particular, they consider a score-based statistic dervied from the marginal means of the observations modelled separately from the correlation structure of the data. The mean of each observation is which is devoid of the of the cluster effect θ ct . The "quasi-score" statistic is then: vector of modified intervention indicators equal to 1 if the intervention was present in cluster c at time t and -1 otherwise, and where n c = t n ct and n ct is the number of individuals in cluster c at time t. G c is a (1 × n c ) vector with elements (∂h −1 ict /∂η ict ) −1 , and V c is an (n c × n c ) covariance matrix for cluster c with non-zero elements off its diagonal. As an example, if we assume the data are normally distributed with mean given by Equation (6), identity link function, variance σ 2 , and θ ct ∼ N (0, τ 2 ), then the diagonal elements of V c are σ 2 + τ 2 and the off-diagonal elements are τ 2 . More complex structures might include temporal decay in correlation, for example.
We use Θ to represent the parameters of the variance-coviarance matrix.
For the permutation test to be valid the 'nuisance' parameters (µ, β, Θ) must be invariant to permutation [Braun and Feng, 2001]. This means we cannot re-estimate them for each new permutation. In practice the maximum likelihood estimates of these parameters are used to construct the test statistic, so that we use the estimates: for the linear predictor under the null H 0 : δ = δ * . Estimating Θ is more difficult, however, particularly when the number of clusters is small [McNeish andStapleton, 2016, Watson et al., 2021]. As an alternative to (11) we can replace G c V −1 c with a (1 × n c ) vector of ones 1 c : so that the sum of residuals is weighted only by the size of each cluster or cluster-time period. Braun and Feng [2001] provide proof of the optimality of the quasi-score (11)  For the purposes of multiple testing we use a studentized version of (13): where the terms on the right-hand side have been evaluated at δ = δ * . In the absence of studentization, the variances of the test statistics are not scale-free and depend on, among other things, the null hypothesis being tested so that different tests will have different power under stepdown testing [Romano and Wolf, 2005a]. The lack of balance is particularly consequential for the construction of confidence sets discussed in the next section. While confidence sets constructed on the basis of stepdown methods will have joint coverage of 1 − α, without balance the individual coverage probabilities of each interval will differ, perhaps substantially [Romano and Wolf, 2005b].
The test statistic described in this section can be used with the stepdown procedure of Section 2.2 to provide a basis for generating p-values for multiple hypothesis tests across multiple outcomes with a nominal FWER. However, a point estimate and p-value is often not considered sufficient for providing inferences from trials. The p-value can provide relatively little indication of the magnitude of sampling variation associated with a particular parameter, and often confidence intervals are expected or required.

Confidence sets and multiple testing
The multiple testing problem extends to the construction of simultaneous confidence intervals or a "confidence set". Let the parameters of interest be δ j with associated confidence to the FWER, we want appropriate control of the coverage of the 100(1 − α)% confidence set such that the process produces confidence sets with the property: we refer to this as 'family-wise coverage'. If we construct 100(1 − α)% confidence intervals independently then the probability that at least one interval in the set excludes the true value can significantly exceed α. There have been some attempts to construct exact confidence sets for parameters based on the stepdown procedure [Romano and Wolf, 2005b].
For example, Guilbaud [2008], extending the proposal of Hayter and Hsu [1994], uses the acceptance/rejection of null hypotheses by the stepdown procedure as a basis of determining upper or lower limits of confidence intervals if we conclude they are strictly negative or postitive, respectively. However, these procedures can only provide information on the upper or lower bound respectively -the other end of the interval is infinity -so they provide little extra information on the extent of sampling variation beyond the p-value.
As an alternative, consider for a moment, a single parameter δ 1 . Its 100(1 − α)% confidence interval is [L 1 , U 1 ]: for any value δ * 1 inside this interval the null hypothesis H 1 : δ 1 = δ * 1 will not be rejected in favour of the two-sided alternative H 1 : δ 1 = δ * 1 at the α level. The question is then how to find the values of L 1 and U 1 efficiently. As suggested by Braun and Feng [2001], one could iteratively perform a series of permutation tests to identify the limits as U 1 = sup{δ * 1 : do not reject δ 1 = δ * 1 } and L 1 = inf{δ * 1 : do not reject δ 1 = δ * 1 }. However, this procedure is inefficient, particularly when testing multiple parameters: if there are M permutations per test and J outcomes, then for each increment in U we must calculate JM permutation test statistics and perform the stepdown procedure. Moreover, since the test statistic and its randomisation distribution depends on the values of the other null hypotheses being tested, a very large number of combinations of values of the parameters must be tested to ensure we have identified with reasonable certainty the limits of the confidence set. Garthwaite and Buckland [1992] developed a method for searching for confidence interval endpoints efficiently, which Garthwaite [1996] later adapted for use with randomisation tests. Their method is based on the search process devised by Robbins and Monro [1951], who developed a stochastic approximation procedure to find the α-quantile of a particular distribution. Multivariate Robbins-Monro processes follow the same procedures as their univariate equivalents [Ruppert, 1985]. For our multiple testing scenario the upper limits to the confidence set correspond to where all hypotheses H j : θ j = U j for j = 1, ..., J are all rejected in favour of the two-sided alternative with a FWER of α but for any smaller values of U j not all hypotheses are rejected, and equivalently for the lower limits.
Returning to Equation 4 and following the logic of the stepdown procedure, if the FWER rate for this family of hypotheses is preserved then P r(|T K j (aX)| ≥ |T |j| (X)|) = α (and equivalently P r(|T K j (aX)| < |T |j| (X)|) = 1 − α) where T K j (aX) is the largest of the permutational test statistic corresponding to the hypotheses H |j| to H |J| . This provides a probabilistic basis for the search procedure based on single iterations of a randomisation test.
At the qth step of Q steps total, we have estimates of the upper confidence interval limits of our J parameters u q = [u 1q , u 2q , ..., u Jq ]. We generate the set of test statistics T j (X)| δ=u jq , which correspond to the null hypotheses H j : θ j = u jq . We then generate a single permutation of a randomisation test for the same hypotheses |T j (a q X)| δ=u jq . The probability that |T j (a q X)| δ=u jq is less than |T j (X)| δ=u jq is α q . The estimates of the upper limits are then updated as: where s j is the "step length constant". The mean increment of each estimate at each step is α q s j (1 − α)/q − (1 − α q )s j α/q = s j (α q − α)/q so that it is clear that each step reduces the expected distance from U j . Similarly for the lower limits, the updating rule is: The step length constants are s j = k(u jq −θ j ) and s j = k(θ j − l jq ) for the upper and lower limits, respectively, whereθ j is a point estimate of the parameter and: where z α is the α-quantile of the standard normal distribution. Under relatively weak regularity conditions this process converges in probability to the correct confidence limits as l → ∞ Buckland, 1992, Ruppert, 1985].

Computation
An R package to execute the analyses described in this paper is available from CRAN as crctStepdown.

Simulation studies
To examine the FWER and simultaneous coverage of the procedures outlined in this article we conduct a series of simulation studies. We start with the simplest relevant design and move to more complex constructions. A brief description of each simulation is described here, a full description can be found in the Supplementary Material.
Each simulation is based on a hypothetical or real two-arm cluster randomised trial of different designs. All outcomes are simulated and modelled using exponential-family models. In all simulations we set the number of individuals per cluster to 20. We consider both when all treatment effects are zero and when only a subset are and also vary the number of clusters per arm. An important parameter for cluster trial design and analysis is the intraclass correlation coefficient (ICC), which is the ratio of cluster-level to total variance. For non-linear models the ICC is not constant as it depends on the particular realisations of covariates and parameter values. We set cluster-level random effect terms to have variance to give ICCs in the range of 0.01 to 0.10. The treatment effect parameters for each simulation are a vector, δ, with length equal to the number of outcomes.
In each simuation, we compare a small and a large number of clusters per arm, with the exact numbers given below. Standard uncorrected test statistics, like an F -test, generally do not provide reliable inference for trials with fewer than around ten clusters per arm [Leyrat et al., 2018, Watson et al., 2021. While (small-sample) corrected test statistics are available with variable performance, permutation tests should provide exact inference and we should expect to see nominal FWERs and coverage with both small and large numbers of clusters [Li et al., 2017]. Uncorrected test statistics should only provide nominal FWER with a large number of clusters and a single 'correct' null hypothesis.
The four simulation models we compare are: • (1) Two-arm, parallel cRCT with two outcomes and no covariates. There are two outcomes that are Poisson distributed and only measured in a single post-intervention period.
• (2) Two-arm, parallel cRCT, with three outcomes, and covariates. The three outcomes are Poisson, Bernoulli, and Normally-distributed. We include one normally distributed covariate. As before, all outcomes are only measured in a single postintervention period.
• (3) Stepped-wedge cRCT with four outcome measures and covariates. Of the four outcomes, three were Poisson distributed and one normally distributed. We include one normally distributed covariate. For all simulations of this type we include 20 clusters and ten 'steps' (two cluster per sequence). This simulation is based on an ongoing trial of the effect of a remote consulting intervention on consultation rates at primary care clinics.
• (4) Two-arm, parallel cRCT with three outcomes and no covariates. The purpose of this final simulation was to examine the effect of both very large numbers of clusters and a correlation between trial outcomes. The three outcomes were multivariate normally distributed with standard deviation terms all equal to one and the correlation between any two outcomes equal to ρ.

Simulation methods and evaluation
Each set of simulations is run 10,000 times. We then estimate the FWER for p < 0.05, which would have a nominal rate of 5%, and also estimate coverage of 95% confidence sets. We also rerun each set of simulations but use instead uncorrected t-statistics to obtain p-values and confidence sets for the parameters of interest to compare to the corrected procedure.
The t-statistics we use are the default ones returned by lme4 using the standard errors estimated from the correctly specified model (i.e. with cluster and cluster-time random effects, where relevant). 95% confidence intervals for each parameter were constructed in the default way ofδ + / − 1.96 × SE. To conduct the simulations we used the R package crctStepdown described above. Table 1 reports the results from simulation experiments (1) to (3). Under all tested conditions the FWER (for p < 0.05) and coverage (95%) were nominal under the proposed corrected procedure. Running time for a single iteration of the procedure ranged from 20 seconds (for the two outcome, single-period simulation) to two minutes (for the four outcome, stepped-wedge design). Figure 1 shows an example of the confidence interval search method output. Table 2 reports the results of simulation experiment (4). As before, FWER and coverage were nominal in all cases.

Results
We compare our results to a "standard analysis" that uses uncorrected t-tests as the  Table 2 shows FWER of 0.05 for 50 clusters per arm and one true null (and equivalently FWER of around 0.15 for three true nulls). In the simulations where there were a small number of clusters rates were not nominal even when there was only one true null, reflecting previous results [Leyrat et al., 2018, Watson et al., 2021. None of the confidence sets had appropriate coverage under the uncorrected method.

Applied example
To provide a real-world example of the the use of the methods proposed in this article, we re-analyse a cluster randomised trial of a financial incentive to improve workplace health and wellbeing in small and medium sized enterprises (SME) in the United King-  (1)

Outcomes
A single primary outcome was specified in the protocol, which was the question "Does your employer take positive action on health and wellbeing?". However, given the potential lack of insight it might provide into the functioning of the intervention, several secondary outcomes were specified to capture the "causal chain" between intervention and employee health and wellbeing. For each of three separate health categories (mental, musculoskeletal, and lifestyle health) employees were asked: 1. whether the employer provided information in this area; 2. whether the employer had provided activities and services in this area; 3. whether the employee had made a conscious effort to improve in this area; 4. whether the employee had attended any groups or activities in this area at work; 5. whether the employee had attended any groups or activities in this area outside of work.
for a total of 15 outcomes.

Re-analysis
The original analysis of the trial took a Bayesian approach. The Frequentist re-analysis we conduct here is principally for illustrative purposes, and so we only take a subset of the data and simplify some of the outcomes. In particular, we take only the main control arm and the high incentive intervention arm to estimate the effect of the high incentive. We focus on the set of secondary outcomes listed above, which we collapse into five separate outcomes; whether the employer provided information across all three health areas, and then whether there was a positive response for any of the health areas for the remaining outcomes, for a total of five outcomes. All outcomes are modelled using a Bernoulli-logistic regression model, following the notation above, with t = 0 for baseline and t = 1 for post-intervention: Table 3 shows the results of an 'uncorrected' univariate analysis using the standard output from R's lme4 package alongside 'corrected' results using the method proposed in this article. We make several observations. The uncorrected analysis would suggest there is likely good evidence that the intervention improved employer provision of information and activities and services, and increased employee taking part at work. However, this conclusion might contradict our understanding of the causal processes since it would seem contradictory for employees to make more effort but not report making more effort. The results corrected for multiple testing appear to be more consistent in that employers appeared to make more effort but the employees did not take up the new services. The effect of the intervention is also more uncertain than suggested by the uncorrected confidence intervals. In particular, the confidence intervals under the corrected method, which are

Discussion
This article presents the first method to provide Frequentist statistics for cluster randomised trials with multiple outcomes that controls for the FWER and coverage of simultaneous confidence intervals. We build on a range of previous work including: randomisation tests for cluster trials [Gail et al., 1996, Gallis et al., 2018, univariate methods for corrections for multiple testing that use randomisation tests Wolf, 2005a,b, 2016], and procedures for estimating confidence interval limits based on randomisation tests Buckland, 1992, Garthwaite, 1996]. Altogether the proposed method can deal with several issues that are common to cluster randomised trials as it allows for multiple outcomes, it can incorporate other features such as restricted randomisation methods, which are often used in trials with a small number of clusters [Watson et al., 2021, Li et al., 2017, and it provides exact inference when there are a small number of clusters, which can lead to non-nominal error rates of standard test procedures and hence confidence intervals with non-nominal coverage. Several small-sample corrections exist that can provide nominal error rates with a small number of clusters [Watson et al., 2021], however there is no obvious way these would be incorporated efficiently into a multiple testing procedure.
We have used optimal methods for: (i) correction for multiple testing; (ii) choice of test statistic; and (iii) confidence interval limit search. However, while each might be optimal for its chosen problem, in combination they may still be inefficient. Since there are no current alternatives for the problem we consider here though, we cannot compare power under different approaches and we did not consider power. These methods are useful for the analysis of cluster trials with multiple outcomes and the treatment effect parameters from the linear predictors of multiple univariate models, however, it is not clear how or if they could be applied to cluster trials with multiple arms. In multi-arm trials there may be one or more outcomes, but clusters may receive different 'doses' or variants of the treatment. There are a variety of treatment effects and null hypotheses of interest including pairwise comparisons between arms and a global joint null, which can be estimated from a single univariate model with indicators for each arm [Watson et al., 2021]. Pairwise null hypotheses in these models do not make statements about the value of the treatment effects in arms outside the pair under comparison as it is left unspecified, so it is not obvious then how a randomisation test could be conducted for the pairwise comparison that is invariant to randomised allocation. The multiple treatment effects of interest in a multi-arm study clearly fall in the realm of multiple testing. Nevertheless, we believe the methods proposed in this article will be a useful tool for the analysis of cluster randomised trials in many cases.

A Simulation details
A.1 (1) Two-arm, parallel cRCT with two outcomes and no covariates For the first simulation we consider a parallel cluster trial with two outcomes measured once in the post-intervention period. Both outcomes are specified as Poisson distributed variables: for j = 1, 2, where θ j,c ∼ N (0, σ 2 θ ). We set µ j = 1 and consider both δ = (0, 0) and δ = (0, 0.5) to compare the FWER under different combinations of true null hypotheses.
Analyses of cluster randomised trials usually consider the intraclass correlation coefficient (ICC) ρ j = V ar(θ j,c )/V ar(Y j,ic ), however the ICC for non-linear models depends on the realised values of the covariates and the parameter values and so will differ between simulations. We choose σ 2 θ to give a reasonable range of ICCs; for this simulation σ 2 θ = 0.2 2 gives ICCs in the range of approximately 0.01 to 0.10.

A.3 (3) Stepped-wedge cRCT with four outcome measures and covariates
For the final set of simulations we reproduce a simulation used in the design of an ongoing stepped-wedge cluster randomised trial. The trial evaluates an intervention designed to enable and facilitate the use of remote (via telephone or smartphone) primary care consulting in a number of low and middle income countries. The intervention was rolled out to 20 primary care clinics in a stepped-wedge design over ten monthly steps comprising two clusters per step. Patients with a set of pre-specified chronic conditions were included in the study as they regularly used the clinic.
The researchers were interested in the effect of the intervention on: (i) the rate of remote consulting (number of remote consultations per patient per month); (ii) the rate of face-to-face consulting (number of face-to-face consultations per patient per month) -it was unknown if remote consultation functions as a substitute, complement, or neither for face-to-face consultation; (iii) the prescribing rate of medication (number of prescriptions issued per patient per month), which was seen as the primary safety endpoint; and (iv) the trust the patient had in the consultation process (evaluated using a 'trust score' from a specialised survey instrument). All four outcomes were therefore specified as 'primary outcomes'. We present analysis of the real data from this trial in Section 6.