Weighted Empirical and Euclidean Likelihood Covariate Adjustment


 Covariate adjustment is often used in statistical analysis of randomized experiments to increase efficiency of estimators of treatment effects. In this paper, we study covariate adjustment based on the empirical and Euclidean likelihoods, and propose weighted versions that arise as natural alternatives. The weighted methods incorporate the auxiliary information that the covariates have equal means among the treatment groups due to randomization. We show that the empirical and the Euclidean likelihoods and their weighted versions are first order equivalent to Koch’s nonparametric covariance adjustment. Allowing the weights to be negative, the resulting pseudo Euclidean likelihood is equivalent to Koch’s method, and its weighted version can be viewed as a weighted version of Koch’s method. In a simulation study, we assess the finite sample properties of the proposed methods. The analysis of a clinical trial data set illustrates an application of these methods to a practical situation.


Introduction
Analysis of covariance is used in statistical analysis of randomized experiments to increase efficiency of estimators of treatment effects and to induce equivalence (with respect to the covariates) of the treatment groups generated by randomization (Snedecor and Cochran, 1980). Nonparametric alternatives include the rank analysis of covariance (Quade, 1967), the analysis of covariance based on general rank scores (Puri and Sen, 1971), and the nonparametric covariance adjustment for categorical data (Koch et al., 1982). The presence of covariate information in randomized experiments is similar to the presence of auxiliary information in survey sampling (Fienberg and Tanur, 1987). Auxiliary information, mostly in the form of known finite population means of auxiliary variables, is frequently used in survey analysis to increase the precision of estimators (Deville and Sarndal, 1992).
Nonparametric covariance adjustment strategies are described by Koch et al. (1998) for continuous, dichotomous, and ordinal data in randomized clinical trials with two groups. The authors described extensions of these methods for stratified studies, for multivariate responses, for missing data problems, and for estimation of nonparametric measures of treatment effects. Tangen and Koch (2001) further extended the methodology for randomized clinical trials with multiple groups. The power of non-inferiority and superiority tests based on Koch's estimator were evaluated by Lesaffre et al. (2002) and Lesaffre and Senn (2003) proposed correction factors to better maintain the nominal levels of the tests on small samples.
The empirical likelihood (EL) was introduced by Owen (1988) as a semiparametric method of inference for general parameters defined by statistical functionals. An advantage of the EL is that it defines an empirical likelihood ratio test (ELRT) which has similar asymptotic properties as the likelihood ratio test (LRT) in the parametric case. The confidence regions obtained by inverting the ELRT enjoy desirable properties: have a data-determined shape, respect the range of parameters, are transformation respecting, and do not require estimation of scale or skewness (Hall and La Scala, 1990). The ELRT statistic was shown to be Bartlettcorrectable by DiCiccio et al. (1991). Applications of EL for randomized experiments were investigated by Adimari (1995), Jing (1995), Liu et al. (2008), and Owen (1991). Weighted EL methods have been proposed by Chen and Sitter (1999), Ren (2001), Glenn and Zhao (2006), and Tsao and Wu (2006). A comprehensive review on applications of EL is provided by Owen (2001) and on EL methods for sample surveys by Wu (2002).
In this paper, we propose covariate adjustment methods for randomized clinical trials that incorporate the auxiliary information of equal means of the covariates among the treatment groups due to randomization similarly to Koch's approach. These methods include the EL, the Euclidean likelihood (UL), and our proposed weighted versions that arise as natural alternatives: the weighted empirical likelihood (WEL) and the weighted Euclidean likelihood (WUL). The UL uses the Euclidean distance instead of the Kullback-Leibler divergence in the optimization problem corresponding to the EL. Allowing the weights to be negative in the case of Euclidean distance, we show that the resulting pseudo UL is equivalent to Koch's method. The weights used for the WEL, WUL, and pseudo WUL incorporate the available auxiliary information. The pseudo Euclidean likelihood (weighted or unweighted) methods have the computational advantage of providing closed-form expressions. We prove that the EL, the UL, the pseudo UL, and their weighted versions are first order equivalent and that the pseudo WUL can be viewed as a weighted version of Koch's method. Our approach is similar to the pseudo EL for sample surveys of Chen and Sitter (1999), Wu and Rao (2006), and Zhong and Rao (2000). Specifically, both our WEL approach and pseudo EL use weights on the components of the EL function. However, while the pseudo EL uses survey weights, we use the weights derived from the auxiliary information. We will show an application of our weighted approach for sample surveys which can be viewed as a weighted version of the pseudo EL of Chen and Sitter (1999) and Wu and Rao (2006).
Although previous work has described the use of Koch's nonparametric covariance adjustment method for multivariate data structures, the focus has not been on multivariate test statistics (Saville et al., 2011;Kawaguchi et al., 2011;Saville and Koch, 2013). As such, the consideration of multivariate test statistics in their own right and the construction of related confidence regions are also novel contributions provided by the current work. Alternative approaches to covariate adjustment have been proposed to estimate treatment effects in clinical trials such as: matching, stratification, inverse probability weighting, and using propensity score as a covariate. For a comparison between these methods, see Elze et al. (2017). For a survey on medical papers that use either covariate adjustment, or subgroup analysis, or baseline comparisons in clinical trials reporting as well as statistical ramifications and recommendations for future practice, see Pocock et al. (2002). For some unexpected results about the efficiency of treatment effect estimates adjusting for covariates for binary outcomes in logistic regression, see Robinson and Jewell (1991). For an application of covariate adjustment for cost-effectiveness data, see Willan et al. (2004). For a principled, practically-feasible approach for covariate adjustment that yields the desired gain of efficiency using semiparametric theory, see Tsiatis et al. (2008). Note that these references consider the case of univariate responses while in our paper we consider multivariate responses.
The rest of the article is organized as follows. We first describe the EL, the UL, and the pseudo UL for a multivariate mean with auxiliary information and introduce the weighted methods as natural alternatives. We then extend the methods for the common multivariate mean and for the adjusted treatment effect models. We describe alternative methods of inference for the treatment effect in presence of covariates such as the Koch's method and multivariate analysis of covariance. In a simulation study, we compare, in terms of the root mean squared errors of the estimators and size distortions of the tests, the proposed methods with alternative methods of inference for treatment effects. We conduct a data analysis of a clinical trial data set to illustrate an application of the methods in practice. We summarize the results developed in this paper in a conclusion and the proofs of theoretical results are provided in the appendix.
2 Multivariate mean with auxiliary information 2.1 Empirical/Euclidean likelihood for the mean Let X 1:n = {X 1 , . . . , Xn} be an i.i.d. sample from a distribution P on R q , and let V X = Var(X). The parameter of interest is µ X = E(X), where X ∼ P. The EL function of inference for µ X is defined as: where µ ∈ R q and P n = {p = (p i ) ∈ R n : ∑︀ n i=1 p i = 1 , p i > 0}. Strictly speaking,l E X (µ) is the log-EL function, but for convenience, we refer to it here and elsewhere as the EL function. Since the objective function is strictly concave, then (1) has a unique solution if µ is in the relative interior of the convex hull of X 1:n . Let p E (µ) = (p E i (µ)) ∈ P n be the solution of (1) assuming that it exists. Owen (1988) showed that Owen (1988) showed that, similarly to Wilks's (1938) theorem for the LRT in parametric models, −2l E Inverting the ELRT, a (1 − α) × 100% confidence region for µ X is given by If instead of the (forward) Kullback-Leibler divergence we use the Euclidean distance in (1), we obtain the Euclidean likelihood (UL) function: Note that (3) is a quadratic optimization problem with inequality constraints which does not have closed form solutions (for a graduate level textbook on convex optimization, see Boyd and Vandenberghe, 2004). If we allow the weights be negative in (3) and maximize over Q n , where Q n = {q = (q i ) ∈ R n : ∑︀ n i=1 q i = 1}, we obtain the pseudo UL function: Here and in other places, we use q to either denote a vector q = (q i ) ∈ Q n or the dimension of X; the notation should be obvious from the context. Strictly speaking, the UL was defined by Owen as the maximizer of (3) over q ∈ Q n (Owen, 2001, section 3.15). However, to be consistent with the definition of EL, we define the UL in (3) as a maximization over p ∈ P n and we calll P X (µ) the pseudo UL function. Using the method of Lagrange multipliers, Owen (2001) derived the solution for q i of (4) in terms of µ (see the appendix): (5) in (4), yields: By (5), Pr(q(µ X ) ∈ P n ) → 1, whereq(µ) = (q i (µ)) ∈ Q n . Hence −2l U X (µ X ) d − → χ 2 q , so that, the EL, the UL, and the pseudo UL are first order equivalent. Note that the pseudo UL has a computational advantage compared to the EL and UL since the weights have closed-form expressions. Another advantage of pseudo UL is that if µ is not in the convex hull of the observations, the EL and UL are undefined (since the feasible set is the empty set), while the pseudo UL and the corresponding weights are always defined.

Empirical/Euclidean likelihood for the mean with auxiliary information
Here and in other places, we use p to either denote a probability vector p = (p i ) ∈ P n or the dimension of Y.
The EL function incorporates the auxiliary information that µ X is known and is defined as: Owen (1991) (2). The maximum empirical likelihood estimator (MELE) of µ Y was introduced by Qin and Lawless (1994) and is defined as: By Theorem 1 of Qin and Lawless (1994), where ⪯ denotes the Lowener partial ordering,μ E Y is first order more efficient thanȲ. If, in particular, rank(V XY ) ≥ p, then Σ Y ≺ V Y andμ E Y is strictly more efficient thanȲ. The UL function is obtained using the Euclidean distance in (7): is the solution (in p ∈ P n ) of (3), then using a similar argument as above, Maximization over q ∈ Q n , by (6), yields the pseudo UL function: Simple algebra shows that (see the appendix): Hence, the EL, the UL, and the pseudo UL are first order asymptotically equivalent.

Weighted empirical/Euclidean likelihood
Our weighted empirical likelihood (WEL) is a natural alternative to the EL for the mean with auxiliary information. The method can be readily extended for more general parameters similarly to Owen (1988). Instead of minimizing the divergence between p ∈ P n and the vector of uniform weights n −1 1n ∈ P n , where 1n ∈ R n is the n-vector of ones, the WEL minimizes the divergence between p andp E = (p E i ) ∈ P n while imposing the same constraints as the EL; specifically, If µ is in the relative interior of the convex hull of Y 1:n , the solution for p i of (11) in terms of µ has the following expression (see the appendix): where ν = (µ X , µ) ∈ R p+q andλ wE (µ) ∈ R p+q satisfies: By Jensen's inequality, and the maximum is attained for p i =p E i , 1 ≤ i ≤ n. Henceμ E Y = argmax µ∈R pl wE Y|X (µ|µ X ), and thus, the maximum WEL estimator and the MELE of µ Y are the same. Theorem 1 shows the large sample behavior of the Lagrange multipliers and of the WEL function in local neighborhoods of µ Y , as well as the asymptotic null distribution of the weighted ELRT. and where If we use the Euclidean distance in (11), we obtain the weighted UL (WUL) function: p U i are defined before (8). Maximizing over q ∈ Q n in (18) and usingq i instead ofp U i , whereq i are defined before (9), we obtain the weighted pseudo UL likelihood function given by: The solutionq w (µ) = (q w i (µ)) ∈ Q n (in q ∈ Q n ) of (19) is given by (see the appendix): Hence, the WEL, WUL, and pseudo WUL are first order equivalent.

Application of WEL in sample surveys
Note that our approach is similar to the pseudo EL for sample surveys of Chen and Sitter (1999), Wu and Rao (2006), and Zhong and Rao (2000). Specifically, both our WEL and the pseudo EL use weights on the components of the EL function. However, while the pseudo EL uses the normalized survey weights, we use the weights derived from the auxiliary information. In this section, we show an application of our method to sample surveys, and in this case, our WEL method can be viewed as a weighted version of the pseudo EL of Chen and Sitter (1999) and Wu and Rao (2006). The weighted EL for sample surveys can be derived in two steps. In the first step, we derive the weights that minimize the divergence to the normalized survey weights (given by the normalized inverse survey selection probabilities) using the auxiliary information; these weights contain both information about the survey weights and the auxiliary information. In the second step, we use these weights (instead of the normalized survey weights) when we define the WEL function for survey data.
Specifically, letd i denote the normalized survey weights corresponding to the sample such that ∑︀ n i=1d i = 1 (for more details, see Wu and Rao, 2006). Then, the pseudo EL function for µ X is given by: Letp E (µ) denote the solution (in p ∈ P n ) of (22). Under regularity conditions (see, e.g., Wu and Rao, 2006), −2l E X (µ X ) is asymptotically distributed as a chi-squared variable scaled by the design effect. In the second step, consider the WEL function for µ Y under the auxiliary information that µ X is known: wherep E = (p E i ) ∈ P n andp E =p E (µ X ). Note that our WEL approach in sample surveys with auxiliary information is different from the pseudo EL approach of Wu and Rao (2006), Chen and Sitter (1999), and Zhong and Rao (2000) and can be viewed as an alternative to the pseudo EL. Notice that our WEL approach to sample surveys blends naturally the weights from auxiliary information and survey weights in statistical models with auxiliary information. The asymptotic and finite sample properties of the corresponding WEL methods in sample surveys for the models considered in this paper is a topic of future research.

Empirical/Euclidean likelihood
where p k = (p k,i ) ∈ P n k , k = 1, 2. Using Lagrange multipliers (see the appendix), we obtain the solutions for p k,i of (24) in terms of µ: Assume that n k /n → π k > 0, where n = n 1 + n 2 . By Theorems 1 and 2 of Qin and Lawless (1994), where, with a slight abuse of notation, Σ X = (π 1 V −1 X1 + π 2 V −1 X2 ) −1 . Note thatμ E X andX k are consistent estimators of µ X ; since Σ X ≺ π −1 k V X k for k = 1, 2, thenμ E X is strictly more efficient than bothX 1 andX 2 . If we use the Euclidean distance in (24), we obtain the UL function: Maximizing (27) over q k ∈ Q n k , we obtain the pseudo UL function: The solutionsq k (µ) = (q k,i (µ)) ∈ Q n k of (28) are given by (see the appendix): (29) for q k,i in (28), yields: Letμ P X be the pseudo MULE of µ X . By (30),

Weighted empirical/Euclidean likelihood
The WEL function is obtained by minimizing the divergence between p k ∈ P n k andp E k ∈ P n k in (24), i.e., wherep E k = (p E k,i ) ∈ P n k are given by (26). Using Lagrange multipliers (see the appendix), we obtain the solutions for p k,i of (32) in terms of µ: Similarly as before, the maximum WEL estimator and the MELE of µ X are the same. Theorem 2 shows the large sample behavior of the Lagrange multipliers and of the WEL function in local neighborhoods of µ X , as well as the asymptotic null distribution of the weighted ELRT.
Theorem 2. Let X (k) 1:n k = {X k,1 , . . . , X k,n k } be two independent i.i.d. samples from the distributions P k on R q , with E(X 1 ) = E(X 2 ) = µ X , where X k ∼ P k , k = 1, 2. Let n = n 1 + n 2 and assume that n k /n → π k > 0 and V X k = Var(X k ) are nonsingular. Then, for every K > 0, and where If we use the Euclidean distance in (32), we obtain the WUL function: (27). Maximizing (38) over Q n k , we obtain the pseudo WUL function:l wP X1,X2 (µ) = max whereq k,i =q k,i (μ P X ) andq k,i (µ) are given by (29). The solutionsq w k (µ) = (q w k,i (µ)) ∈ Q n k of (39) are given by (see the appendix):q and thus, the WEL, the WUL, and the pseudo WUL are first order equivalent.

Empirical/Euclidean likelihood
is the response and X k,i ∈ R q is the covariate corresponding to the ith subject in the kth group, 1 ≤ i ≤ n k , k = 1, 2. The parameter of interest is the treatment effect . The randomization implies that the covariates X k,i have equal means between the treatment groups, i.e., E(X 1 ) = E(X 2 ). The randomization also implies that V X1 = V X2 ; for notational purposes, we do not make this distinction here.
The EL function for (µ X , ∆ Y ) is given by: The profile EL function for where, with a slight abuse of notation, ν = (0, ∆) ∈ R p+q . Using Lagrange multipliers (see the appendix), we obtain: whereη E k (∆) andλ E (∆) satisfy (b) and (c) of (76). (24), andp E k are given by (26), then∆ E Y is the MELE of ∆ Y . By Theorems 1 and 2 of Qin and Lawless (1994), Y is strictly more efficient than∆ Y . Note that, we use here the auxiliary information that the covariate means are equal due to randomization. However, randomization also implies that higher order moments of the covariates are equal between the treatment groups. Our results apply similarly if we consider equality of higher order moments by redefining X to include higher order moments of covariates, and thus, without loss of generality, we only consider the auxiliary information of equal means.
If we use the Euclidean distance, we obtain the profile UL functionφ U Z1, Similarly to (45), (47) over q k ∈ Q n k , yields the profile pseudo-UL function: The solutionsq k (∆) = (q k,i (∆)) ∈ Q n k of (48) are given by (see the appendix): Let∆ P Y be the maximum pseudo UL estimator of ∆ Y . By (50), we obtain (see the appendix): whereq k,i =q k,i (μ P X ) andq k,i (µ) are given by (29). Notice that Koch's estimator has the same expression as (51) (see section 4.4). By (49), then Pr(q k (∆ Y ) ∈ P n k ) → 1, so that

Weighted empirical/Euclidean likelihood
The profile WEL function is given by: wherep E k = (p E k,i ) ∈ P n k satisfy (26). Solving (52), yields (see the appendix): whereη wE k (∆) andλ wE (∆) ∈ R p+q satisfy (b) and (c) of (78). Similarly as before, the maximum WEL estimator and the MELE of ∆ Y are the same. Theorem 3 shows the large sample behavior of the Lagrange multipliers and of the WEL function in local neighborhoods of ∆ Y , and the asymptotic null distribution of the weighted ELRT.
. Assume that n k /n → π k > 0, where n = n 1 + n 2 . Then for every and where The profile WUL function is given by: (27), andμ U X is the MULE of µ X . Maximizing (58) over q k ∈ Q n k , yields the profile pseudo WUL function: whereq k,i =q k,i (μ P X ),q k,i (µ) are given by (29), andμ P X is given by (31). The solutionsq w k (∆) = (q w k,i (∆)) ∈ Q n k of (59) in terms of ∆ are given by (see the appendix): , andq k,i (µ) are given by (29). Substitutingq w k,i (∆) for q k,i in (59), yields the following closed-form expression for the profile pseudo WUL function (see the appendix): Since the pseudo UL is equivalent to Koch's nonparametric covariance adjustment, the pseudo WUL can be viewed as a weighted version of Koch's method.

Alternative estimators for the treatment effect 4.4.1 Unadjusted treatment effect estimation
If the covariate information is not used, the statistic usually employed to estimate ∆ Y is the (unadjusted) treatment effect estimate∆ Y =Ȳ 1 −Ȳ 2 . The Wald test statistic for testing H 0 :

Multivariate analysis of covariance
Consider the multivariate analysis of covariance (MANCOVA) model to estimate ∆ Y : where Y k,i ∈ R p andX k,i ∈ R q denote the response and the (centered) covariate corresponding to the ith subject in the kth group, µ k ∈ R p is the kth group mean, β k,j ∈ R q is the linear effect of the covariates on the jth component of the response in the kth group, B k = (β k,j ) is the p × q matrix of regression coefficients corresponding to the kth group, n k is the number of subjects in the kth group, and ϵ k,i ∼ i.i.d. (0, Vϵ k ) are the random errors, with 1 ≤ i ≤ n k , 1 ≤ j ≤ p, and k = 1, 2. The covariatesX k,i are centered at 0 by subtracting the (overall) mean, i.e.,X k,i = X k,i −X, whereX = n −1 ∑︀ 2 k=1 ∑︀ n k i=1 X k,i and n = n 1 + n 2 . In this parametrization of the model, µ k is the kth group mean response obtained when the covariates are set equal to the (overall) population mean. The parameter of interest is the treatment effect ∆ Y = µ 1 − µ 2 . We rewrite (64) in an equivalent form is the vectorization operator, and ⊗ is the Kronecker product. The least squares estimates (LSE) of µ k and β k are given by (see appendix): whereX k =X k −X. The LSE of ∆ Y is thus given by: SinceX 1 =X 1 −X = (n 2 /n)∆ X andX 2 =X 2 −X = −(n 1 /n)∆ X , then Under regularity conditions Hence∆ M Y is first order equivalent to Koch's estimator.

Simulation models
The computations were done in the R language (R Core Team, 2020). The EL weights were calculated using the R package emplik (Zhou and Yang, 2018) and the UL weights were calculated using the R package quadprog (Turlach and Weingessel, 2019). The WEL method for one sample was implemented using the R package nleqslv (Hasselman, 2019). We implemented the methods in an R library, called WELCAM (Weighted empirical and Euclidean likelihood covariate adjustment methods), which is available from the authors upon request.
In the simulation study, we compared the EL and the UL with Koch's method, the unadjusted method, and MANCOVA in terms of root mean squared errors (RMSE) of the estimators and size distortions of the Wald tests, the ELRT, the ULRT and their weighted versions. We simulated data according to a randomized experiment with two groups, a bivariate response, and a bivariate covariate which may be regarded as a multivariate version of the simulation model of Lesaffre and Senn (2003). Specifically, we generated samples Z (k) 1:n k of i.i.d. draws from a multivariate normal distribution with mean 0 and covariance matrix V Z , where k = 1, 2, p = q = 2, n 1 , n 2 ∈ {50, 100, 200}, and . We also generated data from a centered and scaled Gamma distribution, but since the results were qualitatively similar to the normal case, we include the results for the normal case only; the results for the gamma case are available from the authors upon request.
The parameter of interest is the treatment effect ∆ Y = 0 ∈ R 2 . For each simulation model, we generated S = 10, 000 samples to approximate the RMSE of the estimators and the sizes of the tests.
To improve the small sample performance of the tests, we use an F-distribution approximation of the null distribution of the test statistics. Specifically, letT be a test statistic for the treatment effect with a limiting χ 2 p distribution under H 0 , with two groups of sizes n 1 and n 2 , and a covariate vector X ∈ R q . First, compute the scaled test statistic p −1T ; then, instead of using the critical value of a scaled χ 2 p distribution (i.e., the distribution of p −1 Z 2 , where Z 2 ∼ χ 2 p ), which, with a slight abuse of notation, is the F-distribution Fp,∞, we use the critical values of Fp,n 0 , where n 0 = min(n 1 , n 2 ) − p − q. That is, reject H 0 if p −1T ≥ F p,n0;1−α , where F p,n0;1−α is the (1 − α) quantile of Fp,n 0 . We have also used a small sample adjustment based on a Wishart approximation for the distribution of the sum of independent sample covariance matrices (Nel and van der Merwe, 1986). Since both methods have similar small sample behavior, we did not include the simulation results for the latter. Table 1 shows the RMSE of the unadjusted estimator (∆ Y ), Koch's estimator

Simulation results
. We did not include the results for the pseudo MULE since it is equivalent to Koch's estimator (see section 4.1). Note that∆ Y has the largest RMSE for all sample sizes and correlation coefficients, almost identical values (with slightly larger values for∆ U Y when n 1 = n 2 = 50). While the RMSE of∆ Y are similar for the same sample size and different correlation coefficients, the RMSE of and∆ U Y are smaller for higher correlation coefficients. This shows, as expected, that the adjusted estimators are more efficient for stronger correlations between the response and the covariates.   Table 2 shows the simulation estimates of the sizes of the Wald tests based on the unadjusted estimator (T Y ), Koch's estimator (T K ), and MANCOVA estimator (T M ), the empirical likelihood ratio test (φ E Z ), the Euclidean likelihood ratio test (φ U Z ), the pseudo Euclidean likelihood ratio test (φ P Z ), the weighted empirical likelihood test (φ wE Z ), the weighted Euclidean likelihood ratio test (φ wU Z ), and the weighted pseudo Euclidean likelihood test (φ wP Z ). The critical values of the tests were calculated using the small sample adjustment based on the F-distribution. Note that that all tests maintain their nominal levels very well for all sample sizes, with most empirical sizes between 0.045 and 0.055. The results for the Wald test based on Koch's estimator are slightly different from the pseudo ULRT because we used the unbiased estimators of covariance matrices for the former. The results using the chi-square critical values were inferior for smaller and moderate sample sizes (n 1 , n 2 ≤ 150), while they behaved similarly for larger sample sizes (results available from the authors upon request).

Data analysis
The data set used for this analysis is from the IMPROvE study (Amris et al., 2014), a randomized, controlled, single-center trial with two groups that evaluated the functional and psychological outcomes of a two-week, group-based, multimodal, interdisciplinary rehabilitation program, tailored for patients with chronic widespread pain (CWP). Participants were recruited from the outpatient clinic of the Department of Rheumatology, Frederiksberg Hospital, Denmark. Female patients aged 18 years and older, with CWP diagnosed according to the American College of Rheumatology 1990 definition of widespread pain Wolfe (1990), were eligible for the trial. Exclusion criteria were concurrent psychiatric disorders not related to the pain and other uncontrolled rheumatic or medical diseases that may cause CWP. Based on an a priori sample size calculation, the trial enrolled 192 female patients with CWP randomized 1:1 to either treatment or a waiting list control group, stratified by the baseline level of functional ability as measured with the Assessment of Motor and Process Skills (AMPS) Fisher (2010). One patient withdrew before receiving the result of the randomization, and thus the study sample had 191 female patients diagnosed with CWP who were randomly assigned to either the intervention group (n 1 = 96) or the control group (n 2 = 95). The AMPS is an individualized, observation-based evaluation of functional ability performed by occupational therapists, developed to establish the extent of an individual's ability to execute and complete familiar and relevant activities of daily living (ADL). Two ADL ability measures are provided by the AMPS, one for ADL motor ability (ADLm) and another for ADL process ability (ADLp). The ADLm ability measure is an indication of how much effort and/or fatigue the person demonstrated when performing ADL tasks, while the ADLp ability measure is an indication of how timely and well-organized the person was observed to be during task performance. The AMPS is based on a Rasch measurement model and it provides equal-interval linear measures of ADL ability expressed in logits. The ADLm and ADLp measures evaluated at the 6-months follow-up represent the bivariate response of interest, while ADLm and ADLp measures evaluated at baseline represent the two covariates. The average ADLm at baseline were not significantly different between the treated and controls (1.14 vs 1.13, p = 0.89), and the average ADLp at baseline were also not significantly different between the groups (1.06 vs 1.05, p = 0.96). However, the average ADLm at the 6-months follow-up was significantly higher in the treated group (1.36 vs 1.16, p < 0.01), and the average ADLp at 6-months follow-up was also significantly higher in the treated group (1.13 vs 0.93, p < 0.01).
We performed the data analysis using the covariate adjustment methods described in this paper. The unadjusted treatment effect estimates for ADLm and ADLp were 0.21 and 0.20, respectively, with a Wald test statistic for the overall treatment effect ofT Y = 23.98 (p < 0.01 for WUL vs UL) and we could speculate that the weighted test statistics may have more power. However, a local asymptotic power calculation for the weighted methods and comparison with unweighted methods (Qin and Lawless, 1994;Peng and Schick, 2018) is a topic of future research and beyond the scope of this paper. Figure 1 shows the 95% confidence regions for the treatment effects for the ADLm and ADLp scales obtained by inverting the ELRT, the ULRT, the weighted ELRT (WELRT), the weighted ULRT (WULRT), and the Wald tests based on∆ K Y and∆ Y . The confidence regions are all similar with the exception of inverting the Wald test using the unadjusted estimator which was the largest. This shows that, in agreement with our theoretical and simulation results, the unadjusted method is the least efficient and the empirical and Euclidean likelihood methods and their weighted version behave similarly. Lastly, although the confidence regions produced by different covariate adjustment methods are not the same, they all provide evidence of a significant positive effect of the multicomponent therapy on the ADLm and ADLp scales.

Conclusion
In this article, we proposed weighted empirical and Euclidean likelihood methods of inference for the mean with auxiliary information. In the particular case of randomized clinical trials, the auxiliary information consists of that the covariates have equal means among the treatment groups due to randomization. We studied the asymptotic properties of the proposed methods and showed that they are first order equivalent to their unweighted versions. We showed that allowing the weights to be negative for the UL, the resulting pseudo UL is equivalent to Koch's method. We further showed that the weighted pseudo UL can be viewed as a weighted version of Koch's method. We showed that our method can be applied for inference in sample surveys with auxiliary information and that it blends naturally the survey weights with the weights from the auxiliary information. In a simulation study, we compared the finite sample properties of these methods with alternative methods, such as the Koch's method and multivariate analysis of covariance. We illustrated an application of these methods to the practical situation of data analysis of a clinical trial data set.
Based on the results of our simulation study and its clear computational advantage, we recommend the use of Koch's method of inference for the treatment effect in randomized clinical trials. We have investigated additional simulation scenarios, e.g., when the multivariate response in one group has a normal distribution and in the other group has a scaled gamma distribution, and Koch's method has performed the best overall (results available from the authors upon request). Future research needs to evaluate the performance of these methods for generalized empirical likelihood (Baggerly, 1998), misspecification of auxiliary information (Giurcanu andPresnell, 2018, 2019), for high dimensional case (Giurcanu, 2016), and for more complex models such as causal inference in presence of confounders as indicated by a directed acyclic graph (Rubin, 1974).

Conflict of interest:
Prof. George Luta is a member of the Editorial Board of Open Statistics. He was not involved in the review process of this article.