CLME: An R Package for Linear Mixed Effects Models under Inequality Constraints

In many applications researchers are typically interested in testing for inequality constraints in the context of linear fixed effects and mixed effects models. Although there exists a large body of literature for performing statistical inference under inequality constraints, user friendly statistical software implementing such methods is lacking, especially in the context of linear fixed and mixed effects models. In this article we introduce CLME, a package in the R language that can be used for testing a broad collection of inequality constraints. It uses residual bootstrap based methodology which is reasonably robust to non-normality as well as heteroscedasticity. The package is illustrated using two data sets. The package also contains a graphical user interface built using the shiny package.


Introduction
Inequality constraints arise naturally in many applications. For example, to evaluate if a chemical is a toxin, a toxicologist may conduct a dose-response study to determine if the mean response is monotonic in dose. More precisely, suppose θ i , i = 1, . . . , p with p ≥ 2, are the mean responses of a chemical corresponding to p dose groups. In this case the null and alternative hypotheses of interest are H 0 : θ 1 = θ 2 = . . . = θ p , and H a : θ 1 ≤ θ 2 ≤ . . . ≤ θ p , with at least one strict inequality (known as the simple order constraint), respectively. Sometimes, when the doses exceed the maximum tolerated dose (MTD), it may result in a dose-related toxicity and the monotonicity is violated causing down-turn at some (unknown) dose i (Simpson and Margolin 1986). In such cases, researchers are interested in testing for an umbrella alternative H ai : θ 1 ≤ θ 2 . . . ≤ θ i−1 ≤ θ i ≤ θ i+1 ≥ . . . ≥ θ p , with at least one strict inequality.
In a multi-center rat uterotophic assay conducted by the OECD (Organization for Economic Cooperation and Development), researchers were interested in studying the effect of exposure to estrogen like compounds in the uterine weights of pre-pubertal rats. They were interested in testing if the mean uterine weights of animals exposed to estrogen like compounds increased in comparison to the uterine weights of control animals (Kanno, Onyon, Peddada, Ashby, Jacob, and Owens 2003). Thus the alternative hypothesis of interest is H a : θ 1 ≤ θ i , i ≥ 2, with at least one strict inequality, known as the simple tree order. Here θ 1 is the mean of the control group and θ i , i ≥ 2, are the means of the treatment groups.
In cancer trials, it is common for researchers to be interested in evaluating a cocktail of two or more experimental drugs in combination, each tried at low, medium and high doses. In such cases, the typical order restriction of interest is the loop order denoted by {θ control,control ≤ θ control,low ≤ θ control,medium ≤ θ high,high } {θ control,control ≤ θ low,control ≤ θ medium,control ≤ θ high,high }, where θ a,b denotes the mean response corresponding to ath dose of the first treatment and bth dose of the second treatment. The above null and alternative hypotheses can in general be expressed as H 0 : Cθ = c and H a : Cθ ≥ c, respectively, where C is a suitable matrix of zeros, ones and negative ones of appropriate order, θ = (θ 1 , θ 2 , . . . , θ p ) and c is a suitable vector of known scalars, for example a vector of zeros. Some examples of C and c are provided later, and an illustration of some common orders is given in Figure 1.
It is of common interest to perform statistical inference under inequality constraints, such as those described above, in a linear mixed effects model setting, especially in the context of a repeated measures design where a researcher may be interested in detecting trends. However, despite the existence of a large body of literature on constrained inference spanning over five decades and three books on testing for order restrictions (Barlow, Bartholomew, Brenner, and Brunk 1972;Robertson, Wright, and Dykstra 1988;Silvapulle and Sen 2005), it was only recently that researchers developed methods for performing constrained inference in linear mixed effects models (Rosen and Davidov 2012;Davidov and Rosen 2011;Farnan, Ivanova, and Peddada 2014). While Rosen and Davidov (2012) and Davidov and Rosen (2011) developed likelihood ratio based methods, Farnan et al. (2014) developed a residual bootstrap based method that is designed to be robust to non-normality as well as to heteroscedasticity. Furthermore, Farnan et al. (2014)'s methodology allows for modeling categorical as well as continuous covariates.
Surprisingly, not even the popular statistical analysis program SAS (SAS Institute Inc. 2011) has the capability to perform tests under general inequality constraints in a linear fixed effects model, let alone in the context of mixed effects models. As demonstrated in Farnan et al. (2014), statistical methods that are specifically designed for testing inequality constraints are expected to enjoy substantially higher power than the usual omnibus procedures (e.g., ANOVA) which are designed for two-sided alternatives. This observation, together with the fact that there does not exist a general software for performing statistical tests under linear inequality constraints in linear mixed effects models, motivates the current work. In this paper we introduce an R (R Core Team 2016) package, called CLME (for constrainted linear mixed effects; Jelsema and Peddada 2016) based on the distribution-free residual bootstrap methodology developed in Farnan et al. (2014) and available from the Comprehensive R Archive Network (CRAN) at https://CRAN.R-project.org/package=CLME. There are several packages in R which offer constrained fixed effects models, including glmc (Chaudhuri, Handcock, and Rendall 2012) and ic.infer (Grömping 2010), but neither of these appear to offer support for mixed models; the present work fills this void. Furthermore, since the methodology is based Simple Order Figure 1: Illustration of order restrictions. Each circle represents a parameter of interest. Inequalities between two parameters (i.e., circles) are provided by the lines. The vertical axis denotes relative magnitude of connected parameters. No relationship (either <, =, or >) is known among parameters that are not connected. A nodal parameter is a parameter whose order relationship with every other parameter is known a priori or given by the hypothesis that is is being tested. For example, θ 3 is the nodal parameter in the umbrella orders.
on residual bootstrap, CLME does not depend on normality or homogeneity of variances for the residuals or random effects.
The rest of the paper is organized as follows: Section 2 provides a brief description of the constrained inference for linear mixed effects (LME) models presented by Farnan et al. (2014). Section 3 describes the contents of the package CLME along with implementation details. Section 4 provides some illustrative examples using the package, and Section 5 concludes the paper with a summary and some comments on planned developments of CLME.

Definition of the model
denote a linear mixed effects (LME) model where Y is the N × 1 response vector, X 1 is a design matrix of order N × p 1 and θ 1 is the corresponding p 1 × 1 vector of coefficients (often treatment effects). X 2 is an N × p 2 known matrix of covariates, θ 2 is the p 2 × 1 vector of regression coefficients, and U is an N × c matrix of known constants (random effects). For simplicity we write X = (X 1 : X 2 ) and U = U 1 : U 2 : . . . : U cq , where : denotes columnbinding and U i is an N × c i matrix, with q i=1 c i = c. We also denote θ = θ 1 , θ 2 and p = p 1 + p 2 .
The random vector ξ = ξ 1 , ξ 2 , . . . , ξ q is c×1, where each ξ i is a c i ×1 vector corresponding to U i , for i = 1, . . . , q. The elements of ξ are independently distributed with mean 0 and covariance matrix T = diag τ 2 1 I c 1 , τ 2 2 I c 2 , . . . , τ 2 q I cq . The residual term is similarly defined with mean 0 and covariance matrix Σ = diag σ 2 1 I n 1 , σ 2 2 I n 2 , . . . , σ 2 k I n k , where i = 1, . . . , k and k i=1 n i = N. Although the above model description and the methodology implemented in CLME allows for fairly general settings, in many applications one may not require the full available flexibility. For example, in most applications it may be sufficient to assume that T = τ 2 I, instead of the general heteroscedastic structure for T described above.
Let C be an r × p matrix so that Cθ represents the linear combinations which are subject to inequality constraints specified by the alternative hypothesis. Thus the hypotheses of interest are given by: such that at least one of the r inequalities is strict. Farnan et al. (2014) suggested the pool adjacent violators algorithm (PAVA) to implement the order constraints (isotonization). We depart from their methodology in that we use the package isotone (De Leeuw, Hornik, and Mair 2009) for isotonization, and in particular the function activeSet. In some cases using active set methodology leads to the same results as using PAVA; though using active sets is a more general approach and enables easy specification of complex order restrictions. This also enables access to a number of solvers for activeSet, including least squares, least absolute deviation, and others.
CLME is designed to implement two general classes of statistical tests. The likelihood ratio type (LRT) statistic (Rosen and Davidov 2012) is the default setting, but the user may instead choose the Williams' type test statistic (Williams 1971(Williams , 1977Peddada, Prescott, and Conaway 2001;Farnan et al. 2014). In both cases, to keep the methodology robust to non-normality and potential heteroscedasticity, the p values are evaluated using the residual bootstrap methodology developed in Farnan et al. (2014). Thus, although our likelihood ratio type statistic is motivated by the likelihood ratio principle under the normality assumption, it does not use the normal theory based asymptotic distribution for the test statistic. Hence we use the phrase "likelihood ratio type test" rather than "likelihood ratio test", and results from CLME will not always align with those of a direct implementation of Rosen and Davidov (2012). Farnan (2011) and later Farnan et al. (2014) investigated the performance of the residual bootstrap based test using the above defined likelihood ratio type test (LRT) statistic and the following Williams' type statistic (W ) under a wide range of distributions and variance structures. For example, Farnan (2011) generated data from a wide range of non-normal distributions including, gamma, log-normal and mixture of normals and different patterns of variances to study the Type I error and power of these tests. A subset of these results were reported in Farnan et al. (2014). The residual bootstrap based methodology, using the above test statistics, generally performed well in terms of Type I errors and power. The Type I errors were generally close to nominal levels.

Performance of residual bootstrap methodology
The Williams' type statistic is defined in general as: where denotes the Schur product of vectors, i.e., a b = (a 1 b 1 , a 2 b 2 , . . . , a r b r ) ,θ 1 denotes the estimator of θ 1 under the inequality constraint of interest, andθ 1 denotes the unconstrained estimator of θ 1 (e.g., the MLE). For a given order restriction specified by C, the contrast matrix B is derived from the largest hypothesized difference(s); in the simple order this is the difference between θ 1 and θ p 1 . The structure of B is similar to that of C, and is further described in the paragraph titled "Constraints" in Section 3.1.

Contents of CLME
In this section we describe the functions included in CLME and provide some notes on their implementation. We start by describing the main function of the package, clme. Afterwards, we detail some of the secondary functions which users may find useful.

Main function
The main function of CLME is clme. This function implements the order restricted residual bootstrap test described in Farnan et al. (2014). Among the arguments listed below, only the formula and the data set are required for the model to run. A series of flowcharts are provided in Appendix A (Figures 8-10) to guide a user through the specification of the arguments for clme.
formula: a formula expression; the constrained effect(s) must come before any unconstrained effects.
data: data frame containing the variables in the model.
gfix: optional vector of group levels for residual variances. Data should be sorted by this value.
constraints: list containing the constraints.
tsf: function to calculate the test statistic.
tsf.ind: function to calculate the test statistic for individual contrasts.
verbose: logical, prints iteration step. Argument can be a vector of multiple logicals; successive elements are passed to further functions.
levels: list to manually specify labels for constrained coefficients.
ncon: the number of constrained terms in the formula; the first ncon terms are constrained.
Several of the arguments to clme require further explanation.
Formula. The formula should be a 'formula' object following typical specification in R.
That is, it is a two-sided expression using a tilde operator to separate the dependent variable from the independent variables. For example, the formula y~x1 + x2 is interpreted as modeling y as a linear function of x1 and x2. Random effects are specified in a similar manner as in the popular lme4 package (Bates, Mächler, Bolker, and Walker 2015), using parentheses. For example, the formula y~x1 + x2 + (1 | u) is equivalent to the previous formula, but includes u as a random effect.
Note that clme runs a model without the intercept (and removes the intercept if it was requested). This is done so that the parameters are means of interest, rather than the intercept and offsets, thereby simplifying the computations involved.

Constraints.
The argument constraints is a list describing the order restrictions using the following elements: order: text string specifying the type of order. Allowed values are "simple", "umbrella", and "simple.tree".
node: numeric indicating which element of θ 1 is the node.
decreasing: logical indicating decreasing order. For simple orders, a decreasing order implies a downward trend. For umbrella or simple tree orders, a decreasing order implies a decrease from the node. See Figure 1 for an illustration.
A: matrix describing the order restrictions in C.
B: matrix of coefficients defining the Williams' type statistic (only necessary if Williams' type test is desired). This corresponds to B in Equation 3.
The values of A and B use the same format as the argument isomat from the function activeSet in package isotone: Each is a matrix with two columns where the rows define a specific constraints. For example, A and B are shown below for the decreasing umbrella order with p 1 = 5 and a node at θ 1,3 (the third element of θ 1 ).
The alternative hypothesis is H a : θ 1 ≤ θ 2 ≤ θ 3 ≥ θ 4 ≥ θ 5 . The first row of A defines the constraint θ 1 ≤ θ 2 , the second row defines the constraint θ 2 ≤ θ 3 , and so on. The entries of the rows are the coefficient indices, and the parameter indexed in the left column is hypothesized to be less than or equal to that indexed by the right column. The B matrix is similarly structured, but defines only the contrasts for the largest hypothesized difference(s). Under the umbrella order, this will be the node compared to the first and last values; the specific form of the Williams' type statistic from Equation 3 is: Not all of the elements of constraints are necessary. There are three general formats by which to pass the constraints to clme.

Specific defaults:
One may specify only the elements order, node, and decreasing. In this case the program will call an internal function to generate the values of A and B. Allowed values for order are "simple", "umbrella", and "simple.tree"; also, the node may be omitted for simple orders. Each of the three elements may also be vector-valued (e.g., order = c("simple", "umbrella")) to test multiple orders.
Custom constraints: Alternatively, the list of constraints may contain A directly. This is particularly useful for specifying custom order restrictions such as loop orders or block orders. When a custom A is passed, the program will ignore any values of order, node, and decreasing. If the Williams' type test is selected, a custom B is also needed.
Unspecified: Finally, constraints may be left unspecified. In this case the program will search for both simple and umbrella orders with all possible nodes, both increasing and decreasing orders. As with the first case, the program will estimate the order using the maximum test statistic of all the tested orders.
When testing multiple orders the test statistic is taken as the maximum of all the tested orders, and the program will note the order which produced this value as the estimated order. The bootstrap null distribution of the test statistic is constructed from all the order patterns under consideration, not just the estimated order (that is, for each bootstrap sample, the test statistic is computed for all candidate orders, and the maximum is taken). For reproducibility, one may use the seed argument to set the seed for the pseudo-random number generator. The output from any custom tsf should be numeric. Output with length greater than 1 corresponds to multiple global hypotheses being tested. This should not be used for testing all the individual constraints from the A matrix, as these are calculated separately using the tsf.ind argument. If desired, the test statistic function should also specify the names attribute of the test statistic, for example naming the contrast. An example of testing multiple global hypotheses is shown in Section 4.2, a reanalysis of data from the Fibroid Growth Study (Peddada et al. 2008).
Homogeneity of variances: gfix. The model described in Section 2 permits flexibility. In particular, both ξ (if random effects are included) and may be modeled under the assumption of homogeneity or heterogeneity of variances. Currently, each random effect term is modeled with a separate variance component. The argument gfix defines groups for the residual variance(s). By default, the data are modeled with a single residual variance. If gfix is supplied, then each group of gfix is modeled with a separate residual variance. For example if the variable x1 contains treatment groups, then gfix = x1 will produce a residual variance for each treatment group.
The output of clme is a list with elements: • theta: vector of estimates of fixed-effects coefficients, θ.
• theta.null: vector of estimates of θ under the null hypothesis.
• cov.theta: the covariance matrix of the unconstrained estimates of θ.
• ts.glb: test statistic for the global hypothesis.
• ts.ind: vector of test statistics for each of the constraints (each row of A).
• constraints: List containing the constraints (A) and the contrast for the global test (B).
• dframe: data frame containing the variables in the model.
• search.grid: matrix containing the orders to perform a search over.
• cust.const: logical, whether custom constraints were specified, or constraints were generated.
• ncon: the number of constrained effects.
• tsf: function to calculate the test statistic.
• tsf.ind: function to calculate the test statistic for individual contrasts.
• residuals: matrix containing residuals. For mixed models three types of residuals are given.
• random.effects: predicted values of the random effects.
• gfix_group: group membership to define residual variances.
• gran: group sizes for random effect variance components.
• formula: the formula used in the model.
• call: the function call.
• P1: the number of constrained parameters.
• mq.phi: initial values for the random effect variance components.
• order: list describing the specified constraints.
The function clme only fits the model. The user should run the summary method on the fitted 'clme' object to obtain the bootstrap test. The summary method accepts arguments object (an object of class 'clme', the output from function clme) as well as nsim, the number of bootstrap simulations to perform, and for reproducibility, seed, the seed for the random number generator (RNG). Both nsim and seed can be passed to clme, in which case the summary method will use those values. The output of the summary method returns the same fitted object, but appends p values.

Secondary function
Typical use of CLME will involve fitting a 'clme' object, and then calling the summary method for 'clme' objects to conduct inference. In the course of its evaluation, the summary function calls several other functions, of which the primary one is resid_boot. This performs the integral function of obtaining the bootstrap samples for inference. Some of the arguments for resid_boot are equivalent to those of clme, but some additional arguments are provided to allow users greater flexibility. The arguments to resid_boot are: formula: a formula expression, the constrained effect should be the first term on the righthand side.
data: data frame containing the variables in the model. CLME constructs a bootstrap sample as follows: Step 1: Obtainθ 0 , the estimate of θ under the null hypothesis.
Step 2: Compute the observed values of the random effects and residuals. DenoteΨ = UTU +Σ, whereT andΣ are the estimates of T and Σ described in Section 2. Then compute:ˆ Step 3: Standardize the observed values of the random effects and residuals.
where sd(·) denotes the standard deviation of the elements in the vector.
Step 4: Obtain the bootstrap samples. Letν * denote a bootstrap sample ofν and letδ * denote a bootstrap sample ofδ. Then defineˆ * =σ iν * Finally, construct the final bootstrap sample as: In typical use resid_boot will follow the above procedure, and the arguments will be specified appropriately. However, for greater flexibility, the user may submit any numeric vector of the appropriate length for several of the values. In particular: Fixed effects. The vector of fixed-effects coefficients is represented with theta. By default this will beθ 0 , the estimate of θ under the null hypothesis. If the user wishes to center the bootstrap samples at a different location, they may submit an alternative vector for theta. For example, one may specify theta to be the fitted estimate of θ rather than the null estimate (a vector of length p). This will cause the bootstrap samples to be centered at the fitted alternative hypothesis rather than the null hypothesis (these two estimates of θ are available in a fitted 'clme' object, elements theta and theta.null). Or the user may wish to center at an unconstrained point estimate,θ. Note that the only effect specifying theta will have is changing the fixed-effects estimate in Step 4 above. In particular, the observed random effects and residuals will not be affected. Also, null.resids must be set to FALSE for the specified value of theta to be used directly. Otherwise, theta will be projected onto the null hypothesis (so that Cθ = 0) using activeSet from package isotone.
Variance estimates. The argument ssq denotes the vector of estimated residual variance terms, σ 2 1 , . . . ,σ 2 k . If specified, then these values will be used in place ofσ 2 i , i = 1, . . . , k in Step 4 above. If length(ssq) > 1 then the user should specify gfix as the residual variance groups as with the main function clme. The argument tsq is similarly defined, but containing the random-effect variance components τ 2 1 , . . . ,τ 2 q .
Residuals. The observed residuals are denoted by eps. By default these are calculated as shown forˆ in Step 2 above, the residuals from an unconstrained point estimate (e.g., generalized least squares). The user may supply an N ×1 vector containing residuals calculated in some other fashion. Similarly, the observed random effects are denoted with xi, and by default are calculated as shown forξ in Step 2 above. Again, the user may submit a vector with alternative values for the observed random effects.
The only necessary arguments for resid_boot are formula and data, these will be used to obtain the model matrices. If provided, the values of theta, ssq, tsq, eps and xi will be used for bootstrapping. Any of these values that are left unspecified will be computed as described in Steps 1-4 above.

Other package contents
shiny application. The shiny package (RStudio Inc. 2016) offers the ability to develop a graphical user interface (GUI) which implements CLME. A GUI developed in shiny can be run locally or deployed online. This is particularly beneficial to researchers who are not as familiar with R, or programming in general, but wish to use the methods described here. The package CLME includes a shiny application to run clme. After installing the package, a user may run the command shiny_clme() to call the GUI and begin using CLME without any need for further programming. Figure 2 shows the GUI for clme with the arguments filled in. This example uses the rat.blood data set (which was printed to a comma-delimited file). The column headers are: id, time, temp, sex, wbc, rbc, hgb, hct, spun, mcv, mch, mchc, plts, and grp_ord. The final column was added to the data set, it contains the values '0 Hour', '6 Hour', '24 Hour', '48 Hour', and '72 Hour'. This column was added so that the proper order of groups for the constrained effect can be defined. A screenshot of the first 15 rows of this file can be seen in Figure 3.
First the user should browse to the data set (which should be a CSV file with the first row being a header), and then select the desired function arguments. Some arguments are provided by default (e.g., the global test statistic and the type of solver for isotonization).
There are several check boxes to define the order, and then the user must tell the application which variables to use. The variables are identified using either their column number or column letter (e.g., 1 or A). Multiple variables may be separated by a comma (e.g., '1,2,4' or 'A,B,D'), a range of variables may be defined with a dash (e.g., '1-4' or 'A-D'), or a combination of the two can be used. These values should be set to 'None' to indicate no covariates or random effects. Group levels for the constrained effect may not be read into R with the correct order; an extra column may contain the ordered group levels (it may therefore have different length than the rest of the data set). The user is recommended to inspect the bottom of the summary output to verify the ordering of the group levels.   If we wish to model hct, we specify column 8 for the response variable. Then the constrained effect is time, which is in column 2. The other covariates, temp and sex, are in columns 3 and 4, so we type '3-4' (alternatively, we could have typed '3,4', or 'C,D'), and finally the random effect, id, is in column 1. Since this is just for illustration, we specified only 10 bootstrap samples. Figure 2 shows the first panel of output, which contains summary information on the data set: a box plot of the response variable for each level of the constrained effect (if the constrained effect is a factor), and a table of descriptive statistics. The other three panels are: 'Model Summary', which provides the summary printout of the fitted model (several examples are shown in Section 4); 'Model Plot', which provides a plot of the fitted means with indication of significance (e.g., Figure 5); and 'Model Data' prints the data used for the analysis. Figure 4 shows several of the optional parameters in the GUI. By clicking the 'Format Output' checkbox, one can request a confidence interval to show on the plot, adjust the Type I error rate, and force the constrained effect to be a factor (not all file imports will bring in text variables as a factor). The 'Heteroscedasticity' box bring up another input field to select a variable, which allows one to set the gfix parameter of clme. In this case, we have entered '2', which corresponds to the variable time. This has the effect of assuming each time group has a unique variance.
The checkbox 'Define order of constrained groups' will bring up another input field. Here we enter 14, to select the 'grp_ord' variable which was discussed above. Note on the box plot in Generic name Description AIC Computes Akaike information criterion. as Coerces an object to class 'clme'.

BIC
Computes Bayesian information criterion. confint Computes individual confidence intervals for fixed effects parameters. Intervals are centered at the constrained estimates, but use standard errors of the unconstrained estimates. fixef Extracts estimates fixed-effects coefficients, θ (also, fixed.effects and coef). formula Extracts the formula for the model. is Tests if an object is of class 'clme'. logLik Computes the log-likelihood under the assumption of normality.

model.frame
Data frame with the variables in the model. model.matrix The fixed-effects design matrix. nobs Number of observations. plot Produces a plot of the constrained coefficients and denotes statistical significance. print A basic printout of the model results. ranef Extracts predictions of the random effects (also, random.effects). residuals Extract various types of residuals. sigma The residual variance(s). summary Obtains inference (e.g., p values) for objects of class 'clme'.

VarCorr
Estimates of variance components. vcov The variance-covariance matrix of the fixed-effects estimates. Table 1: List of methods currently defined for objects of class 'clme'. All methods extracting some value (e.g., fixef) will also work on the result of summary results of an object of class 'clme'. Figure 2 that the groups are out of order: Without this column, the time periods would be out of order: the '6 Hour' group would be placed between the '48 Hour' and '72 Hour' groups. By specifying this input, we correct the order (output not shown). The checkbox 'Select Control Parameters' allows one to set the seed for the RNG, as well as the maximum number of iterations and the convergence threshold for both the EM algorithm and the MINQUE algorithm.
Methods. The function clme outputs an object of the S3 class 'clme'. The methods available for this class are briefly described in Table 1.

Sample implementation
In this section we demonstrate the use of CLME by applying it to two real-world data sets. Some of the analyses mimic those performed in the original papers but in the context of order-restricted inference. Other analyses are intended to exhibit certain features of the package or compare the available options. We emphasize that these analyses are intended not as a scientific reanalysis, but as an illustration. Consequently some modeling choices, the assumption of homogeneity of variances in particular, are not thoroughly investigated. The data analyzed are included in the package as the data sets rat.blood and fibroid.

Hematologic parameters from Sprague-Dawley rats
In a recent study on the effect the amount of time a sample is stored has on various hematological parameters, Cora, King, Betz, Wilson, and Travlos (2012) conducted a time course study using blood samples drawn from Sprague-Dawley rats. Blood samples from 11 female and 11 male rats were kept at either room temperature 21 • C (the control group) or refrigerated at 3 • C for 6, 24, 48 or 72 hours (see Cora et al. 2012 for more details). Although the authors obtained data on a variety of hematological variables in this repeated measure time course study, we shall focus on hematocrit (HCT) and the white blood cell (WBC) count over time. In the case of HCT we shall illustrate some of the options of CLME while testing for simple order with an increasing trend in time. In the case of WBC we test for simple tree order the mean WBC count in the freezer group was at least as high as that of the 0 hour.
First, we load the package and the data.
with at least one strict inequality, here θ i is the mean corresponding to either 0, 6, 24, 48 or 72 hours. In the second case (Case B), we test for a union of umbrella alternatives. If the null hypothesis is rejected then the algorithm selects the pattern that has largest value of test statistic: vs.
H aB : Thus in (B) the order is unspecified but limited to either umbrella or inverted umbrella orders. Note that simple orders (increasing or decreasing) are a special case of umbrella orders, where the peak is the first or last parameter. The peak or the trough of each umbrella is specified using the specification of node. Case (C) is a repeat of case (A), but there we will assume heteroscedasticity of variances between the time groups.
We initially use the default arguments as far as possible. We use the gender of the rat and the storage temperature of the sample as covariates. The R code to test case (A) is provided below along with the results.
To test case (B) we simply need to omit the constraints from the call to clme. The code and results are given below.

Model based on 1000 bootstrap samples
Observe that the alternative hypothesis in (B) is much larger than the alternative hypothesis in (A). Thus, while the conclusions of tests for (A) and (B) are the same, that the parameters satisfy an increasing simple order, the p value associated with (B) is larger because the alternative hypothesis in (B) is larger than the alternative in (A).
Accounting for heteroscedasticity is simple in CLME. For example, suppose we wish to model each of the time points with a different residual variance. To do this we pass the time groups as the argument gfix, as shown below. We will call this case (C).
R> hct3 <-clme(hct~time + temp + sex + (1 | id), data = rat.blood, + gfix = rat.blood$time, constraints = const, + levels = list(2, levels(rat.blood$time))) R> summary(hct3, seed = 42) Linear mixed model subject to order restrictions Formula: hct~time + temp + sex + (1 | id) -1 Model based on 1000 bootstrap samples White blood cell count. Using the white blood cell count data of Cora et al. (2012), we now illustrate our package for testing a simple tree order. Here the nodal parameter is taken to be the population mean corresponding to the 0 hour group. Since box plots of the residuals (not shown) suggested the variances were potentially equal across the groups, we assume homogeneity of variances. For illustration, in this example we use the Williams' type test statistic. The code and results are below, and the coefficients are plotted in Figure 6.
As an alternative to Williams' type test, we repeated the analysis using the LRT (results not provided) and discovered that the LRT did not reject the null hypothesis at the 5% level of significance (p = 0.252). This discrepancy between Williams' type and LRT is not surprising in view of the simulation study reported in Farnan et al. (2014), which indicated that Williams' type test can be more powerful than LRT in some cases. Peddada et al. (2008) investigated growth rate of uterine leiomyomata (fibroids) in black and white women. Since fibroids are hormonally mediated and there is a drop in estrogen levels as women age, it may be reasonable to hypothesize a reduction in fibroid growth rates. Interestingly, Peddada et al. (2008) reported that for white women the rate of growth of fibroids decreased with age (i.e., simple order with decreasing pattern), whereas they did not find any reduction in the average growth rate of fibroids with age for black women. They defined the three age groups as follows: Young (< 35), Middle (35-44), and Old (≥ 45). We shall now re-analyze their data using the methodology available in our package CLME where the alternative hypothesis for women of each race group is a decreasing simple order. Note that for confidentiality, we use a subset of the data from the Fibroid Growth Study, excluding cases which may be personally identifiable, particularly those with only one fibroid analyzed in the study. This subset of the data represents 240 fibroids on 54 women. The original analysis in Peddada et al. (2008) represented 262 fibroids on 72 women.

Fibroid growth rates
The interest in this case is to test for a simple order for each race using a linear mixed effects model. This analysis serves as a useful illustration of customizing the order restrictions, because it cannot be performed with the default settings of CLME. First we load the data and perform some manipulations to get a factor that we can use. We define the variable race.age to encode the interaction of the race and age variables; with six levels ordered as: young black, middle-age black, older black, young white, middle-age white, and older white.
R> fibroid$initVol <-cut (fibroid$vol, c(0, 14000, 65000, Inf), + c("small", "medium", "large"), right = FALSE) For the interaction between the age and race terms, we require constraints which define a decreasing simple order over age for both blacks and whites, but do not impose any order restriction between blacks and whites. We do this as follows: R> const <-list(A = cbind(2:6, 1:5) [-3, ], B = rbind(c(3, 1), c(6, 4))) R> const To understand the construction of these matrices, recall the parameter vector θ 1 is ordered as: young black, middle-age black, older black, young white, middle-age white, and older white. The groups of the constrained effect are transformed into column indicators, meaning there will be six parameters: where YB denotes 'young black', MB denotes 'middle-aged black', and so on. Hence, the first three elements of θ correspond to the blacks, and the last three elements correspond to the whites.
The A matrix must define the proper order restriction on these elements. The first row defines the constraint θ M B ≤ θ Y B , the first row defines the constraint θ OB ≤ θ M B . The second two rows define similar constraints for the white women. None of the rows define a restriction between any of the first three elements (blacks) and any of the last three elements (whites); hence there is no order restriction imposed between the two races.
To test for a decreasing simple order for both blacks and whites, we must also define a function to compute the Williams' type test statistic of Farnan et al. (2014) for both blacks and whites separately. While the matrix of contrasts is provided above, by default the Williams' type test will take the maximum and report a single test statistic. We require a test statistic for each of these contrasts. This is similar to the function w.stat.ind which calculates the test statistics for the individual constraints. However, submitting tsf = w.stat.ind will test all of the constraints in the matrix constA instead of the contrasts in constB. To correct this, we make a small modification to w.stat.ind. All we have done is replace calls to A with calls to B; this will accomplish our goal of producing a global test for both blacks and whites individually.
We are then ready to run the analysis. For simplicity, we assume homogeneity of variances. Results of the analysis are shown in Figure 7. The code for this figure is given below, since it cannot be produced through CLME.

Comparison with other packages
In this section we provide some comparisons of CLME to lme4, which is a popular package for the analysis of linear mixed models. We will use the rat.blood data, and in particular the response variable hct. We chose this variable because the observed means at each time point happen to be ordered according to a simple order. Hence, the fitted means will not be changed by clme, so the results between clme and lmer will be more comparable.

Summary
In this paper we have introduced the R package CLME for performing statistical tests under linear inequality constraints. It allows the user to choose either the likelihood ratio type statistic or Williams' type statistic. Since it is based on the residual bootstrap methodology it is not dependent on any normality assumption. As demonstrated in the paper, the package is simple to implement with default settings (Section 4.1), and more complex hypotheses (Section 4.2) can be accommodated with relatively little effort.
Due to the flexibility and distribution-free nature of the model, as well as the ease of use, we anticipate that many researchers may benefit from using the order-restricted model implemented in CLME instead of standard ANOVA models. Other than this package, there does not appear to be any software which offers constrained inference for linear mixed effects models.
While the current release is stable, the authors have an interest in further developing the functionality of CLME. There are many potential improvements that we foresee. On the methodological side these include: adding more models, such as logistic models; implementing an automated choice of the number of bootstrap samples (see Jiang and Salzman 2012); allowing for correlated random effects; and adding the ability to perform power or sample size calculations. Furthermore, the software does not currently allow for complex covariance structures for the variance components, such as the AR(1) process, although it may be extended to accommodate such structures. Other projected developments include enabling the program to take advantage of parallel processing to speed up the repetitive calculations for each bootstrap sample. Finally, as noted, shiny offers the ability to create apps, making complex models easily available to researchers without the need to write R code. The included app can be run locally, but shiny apps can be hosted on a server and deployed online. A well-designed and web-based application could put the power and flexibility of CLME at a researcher's fingertips. Future developments include improving the app and deploying it online.

A. Flowcharts to determine arguments
Define data frame and formula.
Determine constraints.
Define test statistic functions: tsf and tsf.ind.
Assume homogeneity of variances for residuals?
Leave gfix blank (defaults to NULL.
Set gfix to be factor containing group membership.
Return to main workflow. Figure 8: Main flowchart to determine arguments (left) and flowchart to determine groups for residual variance (right).

Determine constraints.
Set constraints=list() and specify as elements A and Anull. If needed, also specify B.
All should be matrices.
Leave constraints blank, all default orders tested.
Test one (or more) of the default orders? Simple , umbrella , or simple tree?
To test multiple default orders, each element may be a vector of corresponding type. All combinations will be tested.
Any element may be left blank.
All default values will be tested for that element.
Return to main workflow.
Anything known about the constraints?
Set constraints=list() and specify the elements: order = text string node = numeric scaler decreasing = logical   Figure 10: Flowcharts to determine arguments defining the test statistic.