{'fraction_non_alphanumeric': 0.043793086346277835, 'fraction_numerical': 0.01765610276248574, 'mean_word_length': 4.314575645756458, 'pattern_counts': {'":': 1, '<': 14, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Active Set Algorithm to Estimate Parameters in Generalized Linear Models with Ordered Predictors

In biomedical studies, researchers are often interested in assessing the association between one or more ordinal explanatory variables and an outcome variable, at the same time adjusting for covariates of any type. The outcome variable may be continuous, binary, or represent censored survival times. In the absence of precise knowledge of the response function, using monotonicity constraints on the ordinal variables improves efficiency in estimating parameters, especially when sample sizes are small. An active set algorithm that can efficiently compute such estimators is proposed, and a characterization of the solution is provided. Having an efficient algorithm at hand is especially relevant when applying likelihood ratio tests in restricted generalized linear models, where one needs the value of the likelihood at the restricted maximizer. The algorithm is illustrated on a real life data set from oncology.


Introduction
In many applied problems and especially in biomedical studies, researchers are interested in associating an outcome variable to several explanatory variables, typically via a generalized linear or proportional hazards regression model. Here, the explanatory variables or predictors may be continuous, nominal or ordered. Estimates of regression parameters can be obtained via maximizing a least-squares or (partial) likelihood function. Especially if the number of observations is small to moderate, researchers often encounter noisy estimates of the regression parameters, possibly leading to patterns in the regression estimates that violate the a-priori knowledge of a factor being ordered. In order to improve accuracy of estimates and efficiency of overall tests for associations, it is tempting to use the prior knowledge of orderings in some of the regression coefficients. From a Bayesian perspective, receiving estimators in these type of problems is straightforward using Markov Chain Monte Carlo approaches. Pioneered in a linear model framework by Gelfand et al. (1992), Bayesian approaches have been proposed by Robert and Hwang (1996); Dunson and Neelon (2003); Dunson and Herring (2003). We also refer to the discussion in the latter two papers. To use Gibbs sampling to get the ordered predictor estimate in logistic regression, Holmes and Held (2006) combine the approach in Gelfand et al. (1992) with an auxiliary variable technique. Note that using e.g. flat priors on the regression coefficient vector β it is straightforward to show that the maximum a posteriori estimator is equal to the constrained MLE introduced in Section 2. Although conceptually straightforward, the implementation of these Bayesian approaches is not without fallacies. To not only get point estimates but also assess whether parameters are equal or strictly ordered across level of predictors, one needs to borrow from more frequentist approaches and "isotonize" unconstrained parameter estimates (Dunson and Neelon, 2003). Only then one can accommodate "flat regions", i.e. successive estimates for ordered levels that are equal. Although there exists vast literature on frequentist estimation subject to order restrictions (Robertson et al., 1988), estimation in the specific regression model discussed here has gained surprisingly little attention (Mukerjee and Tu, 1995). This may be due to the fact that setting up algorithms in these type of problems is generally difficult (Dunson and Neelon, 2003), and requires approaches that need to be adapted to specific problems, necessitating a vast literature for numerous cases of order restricted estimation. We mention Dykstra and Robertson (1982); Matthews and Crowther (1998); Jamshidian (2004); Tan et al. (2007); Taylor et al. (2007), or Balabdaoui et al. (2009) discussing computation of order restricted estimates in specific regression problems, and Terlaky and Vial (1998); Balabdaoui and Wellner (2004) or Rufibach (2007) for estimation of probability densities under order restrictions. Additionally, generalizations of the pool-adjacent-violaters algorithm (PAVA) to inclusion of continuous isotonic covariates are discussed in Bacchetti (1989); Morton-Jones et al. (2000); Ghosh (2007); Cheng (2009) in the context of "additive isotonic regression". Estimation in this type of model is usually performed using the cyclical PAVA in connection with backfitting. However, note that we are not in this genuinely semiparametric setting, but rather the number of levels of an ordered factor is given a priori and remains fixed for any number of observations. Recently, a type of algorithm, which has been around in optimization theory for some decades (Fletcher, 1987), has gained considerable attention in the statistical literature: active set algorithms. Dümbgen et al. (2007) use and generalize such an algorithm to compute a log-concave density not only from i.i.d. but even from censored data. An algorithm similar in spirit is the support reduction algorithm discussed in Groeneboom et al. (2008). The latter authors apply it to the estimation of a convex density and to Gaussian deconvolution. A slight generalization of the support reduction algorithm is used to estimate a convex-shaped hazard function in Jankowski and Wellner (2009). Beran and Dümbgen (2009) extend active set algorithms to the estimation of smooth bimonotone functions. They illustrate their algorithm on regression with two ordered covariates, so also treating the example dealt with in this paper. However, Beran and Dümbgen (2009) only consider least squares or least absolute deviation estimation, and at most two ordered factors. In this paper, we propose an algorithm for an arbitrary number of ordered factors, and we also provide a characterization of the solution.
A key feature of an active set algorithm is, that although iterative, it terminates after finitely many steps, and that the solution is finally found via an unconstrained optimization. This implicitly implies that, as opposed to some Bayesian approaches (Dunson and Neelon, 2003), the active set algorithm is not hurt if estimates of subsequent levels turn out to be equal. In Section 2 we show that the estimation of a regression function in generalized linear models (GLM) under the above ordered factor restriction can be easily performed using such an active set algorithm.
Optimal scaling A reviewer drew our attention to optimal scaling, where one seeks to assign numeric values to categorical variables in some optimal way, see e.g. Breiman and Friedman (1985); Gifi (1990); Hastie and Tibshirani (1990), or applied to modeling interactions in Van Rosmalen et al. (2009). In Gifi (1990, Section 2) categories of the original categorical variables are replaced by "category quantifications", and from then on the variables are considered to be quantitative. Note that in the approach discussed in this paper, one does not necessarily look for an optimal transformation, but rather imposes a priori knowledge on a given ordered predictor. In the example analyzed in Section 9 it seems plausible that a higher tumor or nodal stage is associated with a higher risk of experiencing a second primary tumor.
Ordered predictors While the treatment of quantitative and grouped predictors in regression models is straightforward, we briefly review alternative approaches that can be applied to deal with an ordered explanatory variable z. Let us assume the levels of z are coded as 1, . . . , k where k ≥ 2 and the levels are increasingly ordered, i.e. 1 ≤ . . . ≤ k. The most straightforward way to incorporate z as a predictor is simply to ignore the information about the groups and consider it a quantitative variable. This approach implicitly assumes that the group levels represent a true dimension, with intervals measured between adjacent categories that correspond to the chosen coding. If the ordinal values are arbitrarily assigned rather than actually measured, the regression coefficient is then difficult or impossible to interpret. Supposedly the most prevalent approach to incorporate an ordered predictor z in a regression model is to introduce k − 1 dummy variables z 2 , . . . , z k where z i = 1{z = i}, i = 2, . . . k. This approach ignores the additional knowledge of z having ordered levels, entailing that the estimated parameters β 2 , . . . , β k corresponding to the above dummy variables may not be increasingly ordered. This is especially relevant in small sample studies, where noisy estimates may confuse the proper order of dummy variable coefficients. To simplify interpretation of models, especially when interactions are to be incorporated, researchers sometimes resort to dichotomizing a grouped factor, i.e. introducing only one dummy variable z 1 = 1{z ≤ l}, for some 1 ≤ l < k. Here, the additional knowledge about the ordered levels is not used and may cause a substantial loss of predictive information (Steyerberg, 2009, Section 9.1). Another choice may be polynomial contrasts. One then introduces new variables z i = i 2 {z = i}, i = 2, . . . , k. To avoid correlated estimators β i and therefore mutually dependent tests when doing variable selection, researchers generally prefer to modify the design matrix in order to get orthogonal polynomial contrasts. The function as.ordered() in R (R Development Core Team, 2009) does this by default. Gertheiss and Tutz (2008) proposed a ridge-regression related approach to perform regression with ordered factors. Consider the predictor z with ordered categories 1, . . . , k and the linear regression model (1) For simplicity, we do not consider an intercept and only one ordered factor. The vectors z j = (1{z i = j}) n i=1 , j = 2, . . . , k are vectors of dummy variables corresponding to the levels of z, Z is the n × (k − 1) design matrix with the z j 's as columns, y ∈ R n is the response and ε an i.i.d. noise vector where ε i ∼ N (0, σ 2 ). Note that for reasons of identifiability, Gertheiss and Tutz (2008) assume β 1 = 0 and therefore omit β 1 z 1 in (1). Instead of maximizing the original likelihood ℓ(β) over β, Gertheiss and Tutz (2008) instead propose to maximize a penalized version of ℓ: (2) Here, λ > 0 is a tuning parameter. The solution to (2) can be explicitly computed as for a fixed and specified matrix Ω. The idea is that y is assumed to change slowly for adjacent categories, a property of β GT that is "encouraged" by the shrinkage estimator (3). However, note that β GT may still contain two adjacent estimates β i , β i+1 such that β i < β i+1 , a somewhat undesired feature in this setting. Furthermore, if we choose j = 1 as our reference level (and therefore implicitly assume that β 1 = 0), it seems reasonable to demand for the estimated coefficients that they are all positive, what is not ensured by using (3). Finally, further considerations are necessary to determine the tuning parameter λ.
Consider Setting (1) as before. In this paper, we introduce an algorithm to solve the following problem: Maximize ℓ assuming that β 1 = 0 and under the constraint that so that we receive non-negative and adequately ordered estimated parameters β 2 , . . . , β k for the factor levels. This approach is appealing since the available knowledge (or our "prior belief") is precisely exploited. Furthermore, constraining the space of allowed parameters can be interpreted as regularizing the estimator, implying higher accuracy of the constrained estimate (Dunson and Neelon, 2003). This is especially relevant in small samples. As can be seen from (4), we can estimate parameters for an ordered factor such that the constraints 0 ≤ β 2 ≤ . . . ≤ β k are enforced, unlike in Gertheiss and Tutz (2008) where the violation of these inequalities is only penalized. In the latter approach the violation of the first of the above inequalities, the non-negativity constraint, is not even penalized. In addition, our estimator is fully automatic, i.e. no arbitrary choices such as the coding of levels, the determination of a cutoff to pool levels, or the selection of a tuning parameter (such as λ above) or bandwidth are necessary.

Testing in order restricted models
There is a vast literature on likelihood ratio testing in models under linear equality and inequality constraints. For a discussion and further references on (exact) testing under restrictions in the ordinary linear regression model see Perlman (1969), Wolak (1987) and Shapiro (1988). Silvapulle (1994) and Fahrmeir and Klinger (1994) generalize these results to generalized linear models, especially logistic and Cox regression. As can be seen from (14) below, any likelihood ratio test (LRT) is constructed as the difference of the likelihoods at the unrestricted and the restricted maximizer of the (partial) log-likelihood function, which entails that one needs an algorithm to compute the restricted maximizer. Silvapulle (1994, Section 4) describes an ad-hoc approach to find the constrained estimators. However, his algorithm is non-standard and tedious to apply (Silvapulle, 1994, p. 856). The active set algorithm described here is a general framework able to tackle general optimization problems under constraints and therefore able to compute the restricted estimators in the above mentioned tests very efficiently. This facilitates the application of LRTs in this type of problem.

Statistical inference and asymptotics
Typically, deriving asymptotic properties of shape-constrained estimators is hard, but the starting point in all these problems (Groeneboom et al., 2001;Balabdaoui and Wellner, 2004;Dümbgen and Rufibach, 2009) is a characterization of the estimator, since all the estimators are defined as maximizer of some rather involved function. The most prominent example of a theoretical treatment of a shape constrained estimator via its characterization is the greatest convex minorant that characterizes the estimator of a monotone density (Grenander, 1956). In Section 6 we characterize the solution in our problem. Besides being the starting point for a more thorough analysis, a characterization also allows to check whether an algorithm actually delivers the correct solution.
Our contribution We propose an active set algorithm to find estimators in GLMs with ordered predictors. The estimators strictly comply with the constraints and are found very efficiently, and in a finite number of steps. For identifiability reasons, most regression approaches assume that the coefficient corresponding to the lowest level of an ordered factor is equal to 0. Our approach ensures that all coefficients corresponding to higher levels are in fact non-negative as well. In addition, neither the estimator nor the proposed algorithm needs a tuning parameter. Having an efficient algorithm at hand that provides restricted estimates facilitates the application of LRTs to check whether ordered predictors should be included in the model. In addition, we provide a characterization of the estimator. This serves (i) as a benchmark to verify that the algorithm indeed delivers the maximizer, (ii) gives some insight in the structure of the estimator and (iii) marks the starting point for a more thorough (asymptotic) analysis.
Organization of the paper A general formulation of the problem is given in Section 2. Some examples of GLMs that illustrate our new approach are discussed in Section 3. A description of the active set algorithm adapted to our problem is given in Section 4. There exist special cases of the problem that allow one to find the linear regression estimator β 1 more easily than using the active set algorithm, discussed in Section 5. A characterization of the solutions is given in Section 6. Some indications on statistical inference are provided in Section 7. Literature on likelihood-ratio testing to check whether an ordered factor should be included in the model is briefly discussed in Section 8. A real data example from oncology is analyzed in Section 9. Finally, a more technical description of the algorithm and proofs are postponed to the Appendix.

Setup
We consider the general regression problem of modeling an outcome y ∈ R based on some feature vector w ∈ R p . Therefore, we are given a set (y i , (w ij ) p j=1 ) of observations, for i = 1, . . . , n. Write . . , n. The predictors are denoted by w ·j = (w ij ) n i=1 for j = 1, . . . , p. Throughout the exposition, n and p are considered to be fixed. In general, for given y and W, we seek to maximize a real-valued concave criterion function yielding an estimated parameter vector β ∈ R p . Note that to define our estimator and to derive the characterization in Section 6, a model needs no further specification that goes beyond the function L. Ordinary, i.e. unordered, factors are assumed to be already coded as dummy variables, so they are considered quantitative. If an intercept is to be taken into the model, we simply assume it to be a quantitative variable of all 1's. Let c denote the number of quantitative predictors and suppose that the last f predictors w ·j , j = c + 1, . . . , p are ordered factors, each with k j levels (so c = p − f ). Furthermore, the coding is assumed such that w ij ∈ {1, . . . , k j }, i = 1, . . . , n, where a higher number corresponds to a "higher" level of the ordered factor w ·j . Introduce the sets of indices J c,p = {c + 1, . . . , p} and L j = {2, . . . , k j } for j ∈ J c,p . Clearly, the case c = 0 (no quantitative variables in the model) is not excluded. However, we assume to have at least one ordered factor, i.e. f ≥ 1 which immediately implies p ≥ 1. In order to respect the ordinal character of each of the factors w ·j we estimate β based on a new data matrix X ∈ R d . This latter matrix is obtained via modifying the original data matrix W by adding dummy variables for the levels ≥ 2 of the ordered factors. We then constrain optimization of the updated functional L = L(y, X, β) to the constrained space of parameters Here, β j, l is the coefficient of the dummy variable corresponding to the level l of the j-th ordered factor, and k = ((0) c i=1 , k c+1 , . . . , k p ) ∈ R p . For ease of notation, we define B = B(c, p, k). Constraining estimation to B ensures that the estimated parameter corresponding to a "higher" level of an ordered factor is at least as large as those of "lower" levels and all estimated parameters are non-negative. Note that our approach also adds something new if we have an ordered factor with only two levels (note that we always lose the level attributed to the baseline), namely that β j,2 ≥ 0 for this ordered factor.

Examples
We briefly specify the GLMs we provide algorithms for. Extensions to other criterion functions are straightforward.
Linear regression Here, y ∈ R n and we estimate β via maximizing the criterion function ℓ n,1 over all β ∈ B. This latter function is defined as Here, x i· denotes the i-th row vector of X. We emphasize that given L there is no need to further specify a model for the data.
Logistic regression In this case, y ∈ {0, 1} n . Using maximum likelihood estimation (MLE) we obtain the log-likelihood function Clearly, T i are the failure times (possibly unobserved), C i the censoring times, δ i = 1{event has happened} and x i· is the feature vector as before. If we introduce the observed time V i = min{T j , C j } for each unit, let denote the number of individuals at risk after time T i , i = 1, . . . , n. The partial likelihood according to Cox (1972) . . , n and letting t 1 < . . . < t D be the observed (assumed to be distinct, for simplicity) event times, we then easily deduce the log-likelihood function: where α (s) is the above expression belonging to the s-th failure time, s = 1, . . . , D.

Properties of the maximization problems Let us introduce the constrained
and the unconstrained maximizers. The conditions on fixed response y and design matrix X under which η i exist and are unique in logistic regression are well studied (Albert and Anderson, 1984;Santner and Duffy, 1986). Silvapulle and Burridge (1986) specify necessary and sufficient conditions for the MLE to exist in logistic and Cox regression. Since the set B is a closed convex cone, the estimators β i exist and are unique for i = 1, 2, 3 at least under the same conditions as those for η i . Conditions for consistency and asymptotic normality of MLEs in GLMs are provided in Fahrmeir and Kaufmann (1985). In this paper we assume that our design matrix X is such that ℓ n,i is concave and coercive for i = 1, 2, 3.

Active set algorithm to compute β i
In Fletcher (1987) an active set algorithm is described, a useful tool for constrained optimization problems. In connection with likelihood ratio tests (see Section 8) we came across Silvapulle (1994). In Section 4 of this latter paper, it seems as if a version of the active set algorithm is described. However, instead of directly computing the "active set" in each iteration (see below), a crude and computationally expensive "all-subset search" is proposed. In the context of mixture models, the algorithm discussed by Groeneboom et al. (2008) can also be interpreted as a variant of an active set algorithm. In Section 3 of Dümbgen et al. (2007) the general principle of active set algorithms is described in detail, complemented by a discussion of its validity. Here, we therefore limit ourselves to the discussion of the main features and points relevant for the application of the active set algorithm to find the β i 's. We briefly sketch the idea of an active set algorithm, and refer to A for a detailed technical exposition of the algorithm for the problem treated here. Let q denote the number of constraints that compose B, ℓ the function to be maximized and β its maximizer, see Section 3. Define for any index set A ⊆ {1, . . . , q} the linear subspace The function φ maps the indices j and l of the dummy variables forming the ordered factors to the number of constraining inequalities, see (15) in A. The crucial assumption for an active set algorithm is that we have another algorithm available that for any A ⊆ {1, . . . , q} (efficiently) computes provided that V(A) ∩ {β : ℓ(β) > −∞} = 0. Subspaces of the parameter space are considered when violations of the initial constraints appear in the algorithm. In this case, the active set algorithm varies A in a deterministic way, until finally β(A) = β. In order to tailor an active set algorithm to a specific problem, the above maximization on a subspace is crucial. In our regression with ordered covariates setting, we show in A (see Table 4) that three types of subspaces have to be dealt with, depending on the specific violation that occurs.
It is important to realize that by design, the main routine of an active-set algorithm does not need a stopping criterion as e.g. Newton-type algorithms. Once the algorithm has identified the set A that corresponds to the solution β, it performs an unrestricted maximization (here, a stopping criterion may be necessary), which at least in the linear, logistic and Cox regression examples is unproblematic. Verification that a given β is the maximizer can be done by means of Theorem 3.1 in Dümbgen et al. (2007). Additionally, since there are only finitely many subsets of A, the algorithm terminates after finitely many steps.

Special case: An almost explicit solution
To be able to state the following results concisely, let us introduce for every ordered factor j ∈ J c,p the set of indices where the equality constraint β j,l ≥ 0 is active: In this section, we restrict our attention to the case of linear regression with only one ordinal predictor. If in addition Z 1 ( β 1 ) = ∅, that means the constrained estimator has only strictly positive entries anyway, then ℓ n,1 simplifies such that β 1 can be found via solving (7).
Lemma 5.1. If c = 0, f = 1 and Z 1 ( β 1 ) = ∅, the estimator β 1 is where for l ∈ L 1 The proof of this lemma is postponed to B. The solution to (7) can easily be computed using the PAVA (Barlow et al., 1972;Robertson et al., 1988). This latter algorithm performs at most n − 1 iterations until the vector β 1 is found. One of the initial motivations to analyze regression with ordered predictors, and the reason why we included this very specific example, was to see whether this simple and appealing structure can be carried forward to the more general problem of more ordered factors and additional quantitative variables. However, since (i) we were not able to construct a generalized PAVA algorithm that solves our problem and (ii) we are not only interested in the least-squares problem but also treat GLMs, we switched to an active set algorithm.

Characterization of the solution
There are two main purposes of providing a characterization of the estimator β: (i) knowing the structure of the maximizer of ℓ allows one to cross-check the validity of the proposed active-set algorithm and to check whether it has found the correct maximizer of ℓ. (ii) It is well-known that in such constrained estimation problems, the key to deriving asymptotic properties of the estimator such as consistency or rate of convergence is a characterization in terms of directional derivatives, see the discussion in Section 1.
To be able to state the following theorem properly, we introduce the function ψ : (c+ 1)× L c+1 , . . . , p × L p → {1, . . . , d} that maps the original indices (j, l) to the column number of the respective dummy variable in X, or equivalently, to the index i that corresponds to the entry of the vector β ∈ B(c, p, k) that corresponds to β j, l . Specifically, this function is for any j ∈ J c,p and l ∈ L j , By ψ −1 we denote the inverse of this function, i.e. the function that maps the position i of the entry of β to the indices j and l. Now, for each j ∈ J c,p let h j be the vector of distinctive strictly positive values of (β j ) j∈L j for every j ∈ J c,p and any β ∈ B. Using these definitions we split any vector β ∈ B into the following blocks: where u = 1, . . . , |h j | for each j. Here, |.| denotes the dimension of a vector a or the number of elements in a set. Note that |B 1 | + | ∪ j B 2,j | + | ∪ j,u B 3,j,u | = d. Using these blocks, we are now able to formulate the characterization of the solution.
Note that the entries of the gradient at the active constraints β i , i ∈ B 2,j , are not needed to characterize the solution since γ always equals 0 at these positions. Furthermore, the theorem immediately implies for u = 1, . . . , |h j | and j ∈ L j . To illustrate Theorem 6.1, consider the following example: For n = 200 observations we generated a dataset with standard normally distributed errors, three quantitative variables, one (unordered) factor (with three levels) and one ordered factor (with eight levels). The model we stipulated to generate the response y was y i = 2q 1i − 3q 2i + 0q 3i + 0f 1i + f 2i + f 3i + 0o 1i + 0o 2i + 2o 3i + 2o 4i + 2o 5i + 2o 6i + 5o 7i + 5o 8i + ǫ i where q ji ∼ N (1, 2) for j = 1, 2, 3 and i = 1, . . . , n = 200, each level of any factor (whether ordered or unordered) has the same number of observations and these are randomly allocated to the observations. Finally, ǫ i ∼ N (0, 4) for i = 1, . . . , n. The resulting (constrained) linear regression estimates are given in Table 1. Note that for comparison we also added columns for the estimator ρ 1 which is computed similarly to β 1 , but without the positivity restriction β 6,2 ≥ 0. For this estimator, a characterization similar to that in Theorem 6.1 can be given using exactly the same approach.
Let us now illustrate Theorem 6.1. For either quantitative variables or dummy variables corresponding to unordered factors (which in our context are conceptually equivalent), the respective entry of the gradient ∇ β 1 is always 0. As for the ordered factor, for the entries where the positivity constraint is active (i.e. the elements in Z 6 ( β 1 )), the gradient has a value which is not used (and not necessary) for a characterization of β 1 . The sets defined above are for the simulated example: (2.19, 4.65, 4.79).
These sets then yield the following inequalities, according to (10) and (11):

Statistical inference
Having shown how to compute estimators β i for i = 1, 2, 3, the question arises how to perform (frequentist) statistical inference in these models. Deriving consistency, rate of convergence and limiting distributions for estimators similar to β i under standard assumptions is known to be non-trivial. It is therefore not clear how to construct e.g. confidence intervals for our estimated parameters of the ordered factor. By using the characterization given in Section 6, one should be able to derive rates of convergence and even the limiting distribution of β as n → ∞ in a suitably specified model, thereby generalizing the results of Brunk (1970) and Wright (1981) for isotonic regression to our more general setting. This, together with a generalization of the likelihood ratio tests introduced in Section 8 to an arbitrary number of ordered factors, is subject to ongoing research. Note that bootstrap is not without fallacies in these type of models, see Kosorok (2008) and Sen et al. (2009).

Testing for the presence of constraints
There is a vast literature on likelihood ratio testing in models under linear equality and inequality constraints. For a discussion and further references on (exact) testing under restrictions in the ordinary linear regression model see Perlman (1969), Wolak (1987) and Shapiro (1988). Silvapulle (1994) and Fahrmeir and Klinger (1994) generalize these results to generalized linear models, especially logistic and Cox regression. Suppose a researcher wants to test the following hypotheses: Note that the estimator under H 0 can be computed via an unrestricted maximization. It corresponds to a maximization using a modified design matrix X with the columns ψ(c+1, 2), . . . , ψ(c+1, k c+1 ) omitted.
Since under H 0 we need to consider an unrestricted estimator, we have to constrain attention either to (i) only one ordered factor or (ii) a test of inclusion of all ordered factors against their entire exclusion from the model. The potential influence of the additional ordered factor(s) on the response is assessed with H 1 .
In notation similar to Silvapulle (1994), the above hypotheses translate to Following the development in Silvapulle (1994), the likelihood ratio test statistic to test the Hypotheses (13) is defined as The distribution of T LR is a mixture of χ 2 distributions. The weights are in principle fully specified, however, in general hard to compute (Wolak, 1987). As a remedy, one can either use exact Monte Carlo weights (Wolak, 1987) or bounds on the p-value for the above test (Silvapulle, 1994, Proposition 1). As can be seen from (14) any LRT is constructed as the difference of the likelihoods at the unrestricted and the restricted maximizer of the (partial) log-likelihood function, which entails that one needs an algorithm to compute the restricted maximizer. Silvapulle (1994, Section 4) describes an ad-hoc approach to find constrained estimators. However, his algorithm is non-standard and tedious to apply (Silvapulle, 1994, p. 856). The active set algorithm described here is a general framework able to tackle general optimization problems under constraints and able to compute the restricted estimators in the above mentioned tests very efficiently.

A real data example
We illustrate our new algorithm using a data set from oncology, initially analyzed in Taussky et al. (2005). The goal of the study was to assess the impact of treatment-and patient-related factors on the risk of Variable Type Levels (first mentioned = baseline) factor Chemotherapy (CT) yes, CT no Radiotherapy (rt) factor concomitant boost (CB), hyperfractionation (HF) Sex (sex) factor female, male Tumor stage (t) ordered factor 1 < 2 < 3 < 4 Nodal stage (n) ordered factor 1 < 2 < 3 < 4 < 5 < 6 Performance status (ps) ordered factor 1 < 2 <"stage greater than 2" developing a second primary tumor (SPT) of the upper aerodigestive tract within three years after initial therapy, in head-and-neck cancer patients. For a subset of 231 patients that had either been observed at least three years without SPT or experienced an SPT before three years, the endpoint SPT 3 = 1{The patient experienced a SPT at 3 years or before} was defined and modeled using multiple logistic regression. The explanatory variables are described in Table 2.
Researchers assume in general that higher tumor stage, nodal stage, and performance status correspond to a higher risk of experiencing a SPT. It seems therefore appropriate to use our constrained estimator in this setting. In Figure 1, the unconstrained and constrained estimators η 2 and β 2 are displayed (dot and triangle, respectively) as well as profile likelihood confidence intervals for η (α = 0.05). Values of the likelihoods were ℓ n,2 ( η 2 ) = −101.6 and ℓ n,2 ( β 2 ) = −102.2. Estimates for quantitative predictors, i.e. those for age, treatment, radiotherapy and sex turned out to be very similar for η 2 and β 2 . On the other hand, the "prior belief" or assumption of non-negative and increasing estimates for the levels of the ordered factors tumor and nodal stage and performance status was violated by the unconstrained estimator η 2 and "corrected" by β 2 . The original analysis in Taussky et al. (2005) focused on identifying factors that influence the occurrence of SPT. Variables were not taken into account as ordered factors, but were dichotomized. For comparison, we also computed the restricted and unrestricted estimates in this setting, see Table 3. It turns out that parameter estimates and corresponding odds ratios (OR) for the two approaches were similar, except for the nodal status. Note that the effect of tumor stage is reversed, compared to the case where we consider all factor levels (and do not only dichotomize), compare Figure 1.

Extensions
It is straightforward to generalize the set B(c, p, k) to for arbitrary real numbers r j, l . Using such a more general parameter space could be beneficial in connection with finding the minimum effective dose in dose-response models. The dose levels would then take the role of an ordered factor (Wang and Peng, 2007). Our new approach easily allows us to incorporate further predictors of any of the three types described in the introduction to model the response.
Modeling a factor with decreasing levels can be achieved by reversing the coding of the corresponding ordered factor, using the algorithm under the constraint of increasing levels and finally re-reversing the order of the estimates in the vector β. Using this approach, it is straightforward to find the solutions for all combinations of possible orderings of, say, three ordered factors. By computing the value of the criterion function for all these resulting coefficient vectors, one can find the one with the lowest criterion value, an approach related to finding a global maximum in the criterion function described in van der Kooij et al. (2006). Generalizations to further criterion functions, such as other GLMs or least absolute deviation regression with ordered covariates, are straightforward. As for the latter problem, we suggest smoothly approximating the not everywhere differentiable criterion function, as previously discussed in Beran and Dümbgen (2009).

Acknowledgments
The initial motivation for this research grew out of discussion with Lutz Dümbgen while preparing exercises for his lecture "Optimization" during summer semester 2006 at the University of Bern. I thank Leonhard Held for discussions about the Bayesian perspective of the problem, Sarah Haile for proofreading the final version, and my former employer, the Swiss Group for Clinical Cancer Research (SAKK), for the permission to use the data of Taussky et al. (2005). R-functions (R Development Core Team, 2009) to efficiently compute β for linear, logistic, and Cox-Regression are bundled in the R package OrdFacReg (Rufibach, 2009) and available from CRAN.

A Details of the active set algorithm
In this section, we complement the description of the algorithm indicated in Section 2. Recall the sets of indices J c,p = {c + 1, . . . , p} and L j = {2, . . . , k j } for j ∈ J c,p . In order to respect the ordinal character of each of the factors w ·j we introduced in Section 2 the new data matrix X by adding dummy variables for the ordered factors, such that X = w ·1 , . . . , w ·c , x ·ψ(j,l) l∈L j ; j∈Jc,p for dummy variables The function ψ is given in (8). With the above version of coding, l = 1 is considered the reference level for every ordered factor w ·j and the resulting design matrix X is now an element of R n×d where Again, we denote by x i· the i-th row of X, i.e. the values of the "dummyfied" predictors for the ith observation. In order to respect the ordinal character of each of the factors w ·j we then constrain optimization of the updated functional L = L(y, X, β) to the space of parameters B(c, p, k) given in (5). We write ℓ as placeholder for any of the functions ℓ n,1 , ℓ n,2 , or ℓ n,3 (for ease of notation we omit the dependence on n) and the aim is to find for given response vector and matrix of predictors the vector To fit the constrained maximization problem (6) into the framework of Dümbgen et al. (2007), we write the set B given in (5) as for vectors v i ∈ R d . For ease of notation, we have enumerated the constraining inequalities The function φ : (c + 1) × L c+1 , . . . , p × L p → {1, . . . , q} that maps the original indices (j, l) to the "inequality index" i is given by so that the inequalities can be written as for l ∈ L p and j ∈ J c,p . The vectors v i for any i = φ(j, l) ∈ {1, . . . , q} are received via .
Note that all these vectors are linearly independent. Define for any index set A ⊆ {1, . . . , q} the linear subspace and for β ∈ R d the set A of "active constraints": Maximization on subspace The crucial assumption for an active set algorithm is that we have an algorithm available that for any A ⊆ {1, . . . , q} (efficiently) computes provided that V(A)∩{β : ℓ(β) > −∞} = 0, see Section 4. For simplicity and without loss of generality, fix j = c + 1. Then, for a given β the following situations can cause a non-empty set V(A): Note that the situation v s β ⊤ > 0 for any s = 3, . . . , k c+1 can be treated analogously to Case 2 in Table 4.
To compute the unrestricted maximizer β(A) in the three cases given in Table 4, the strategy is to suitably modify the design matrix X. Precisely, we show how to construct new data matrices X i * and a new corresponding function ℓ i * , i = 1, 2, 3 : R d * → R (here, i stands for the corresponding case in Table 4) for a given A * ⊂ {1, . . . , q} in the three cases of Table 4 such that β(A * ) can be immediately derived from It is crucial to realize that the maximization in (16) is unconstrained and the following arguments show that d i * ≤ d in all considered cases. In what follows, we explicitly state the unconstrained maximization problem, assuming that only the case under consideration is present. Apparent combinations of these basic strategies are necessary in case more than one of the three cases described in Table 4 are present.
Case 1 Writing down the maximization problem (16) explicitly, we get where in general M −i is the matrix M with the i-th column omitted and ℓ 1 * (·, Q) is the criterion function corresponding to ℓ, but based on the design matrix Q.
The function −ℓ n,1 can then be written as N l (β l − m l ) 2 + const(y, X).
The minimum of the latter expression under the constraint β 2 ≤ . . . ≤ β k 1 can easily be found using PAVA. 2 Proof of Theorem 6.1 Before coming to the actual proof, we state a necessary lemma.
Lemma B.1. Let a, b ∈ R n be two vectors having the following properties: Then n i=1 a i b i ≤ 0.
The second term disappears due to (12). As for the first term, we invoke Lemma B.1 where ∇ℓ( γ 2 ) takes the role of a and g that of b to finally deduce that (23) is at most 0. 2 Proof of Lemma B.1 First, note that (17) and (18) immediately imply n i=1 a i = 0.
Repeatedly applying this same trick we finally arrive at By means of (18) and (19) the latter expression remains non-positive. 2